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Abstract 

Factorization  is  important  for  both  practical  and  theoretical  reasons.  In  secure  digital  communi¬ 
cation,  security  of  the  commonly  used  RSA  public  key  cryptosystem  depends  on  the  difficulty  of 
factoring  large  integers.  In  number  theory,  factoring  is  of  fundamental  importance. 

This  research  has  analyzed  algorithms  for  integer  factorization  based  on  continued  fractions  and 
binary  quadratic  forms,  focusing  on  runtime  analysis  and  comparison  of  parallel  implementations 
of  the  algorithm.  In  the  process  it  proved  several  valuable  results  about  continued  fractions. 

In  1975,  Daniel  Shanks  used  class  group  infrastructure  to  modify  the  Morrison-Brillhart  algo¬ 
rithm  and  develop  Square  Forms  Factorization,  but  he  never  published  his  work  on  this  algorithm 
or  provided  a  proof  that  it  works.  This  research  began  by  analyzing  Square  Forms  Factorization, 
formalizing  and  proving  the  premises  on  which  the  algorithm  is  based.  First,  this  research  analyzed 
the  connections  between  continued  fractions  and  quadratic  forms,  proving,  among  other  things, 
that  the  square  of  any  ambiguous  cycle  is  the  principal  cycle.  Then,  the  connection  with  ideals  was 
developed,  requiring  a  generalization  to  the  standard  description  and  formulas  for  multiplication  of 
ideals.  Lastly,  the  connection  was  made  with  lattices  and  minima,  allowing  for  a  generalization  of 
the  formulas  relating  composition  with  distance.  These  results  are  fundamental  to  explaining  why 
Square  Forms  Factorization  works. 

This  research  also  analyzed  several  variations,  including  two  different  parallel  implementations, 
one  of  which  was  considered  by  Shanks  and  one  of  which  is  original.  The  results  suggest  that  the 
new  implementation,  which  utilizes  composition  of  quadratic  forms,  is  slower  for  small  numbers  of 
processors,  but  is  more  efficient  asymptotically  as  the  number  of  processors  grows. 

Shanks’  Square  Forms  Factorization,  including  a  concept  he  called  Fast  Return,  has  been  im¬ 
plemented  in  C  and  Magma  and  some  experimental  runtime  analysis  has  been  done.  A  parallel 
version  in  C  has  been  implemented  and  tested  extensively. 


Keywords:  Factorization,  Quadratic  Form,  Daniel  Shanks,  Continued  Fractions,  Lattice,  Ideal 
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Chapter  1 
Introduction 


The  problem  of  distinguishing  prime  numbers  from  composite  numbers  and  of 
resolving  the  latter  into  their  prime  factors  is  known  to  be  one  of  the  most 
important  and  useful  in  arithmetic,  ft  has  engaged  the  industry  and  wisdom  of 
ancient  and  modern  geometers  ...  the  dignity  of  the  science  itself  seems  to  require 
that  every  possible  means  be  explored  for  the  solution  of  a  problem  so  elegant 
and  so  celebrated. 


C.  F.  Gauss  [8] 

Factorization  is  important  for  both  practical  and  theoretical  reasons.  In  secure  digital 
communication,  the  RSA  public  key  cryptosystem  is  often  used.  The  security  of  this  cryp¬ 
tosystem  depends  on  the  difficulty  of  factoring  large  integers.  In  number  theory,  factoring  is 
of  fundamental  importance. 

The  algorithms  for  factoring  larger  and  larger  integers  quickly  have  developed  significantly 
over  the  years.  There  are  many  ways  of  doing  this,  ranging  from  trial  division  to  the  number 
field  sieve.  This  research  focuses  on  a  specific  underdeveloped  algorithm  introduced  by  Daniel 
Shank  1975,  Square  Forms  Factorization,  but  also  analyzes  two  new  variations,  one  due  to 
a  conjecture  by  Pomerance  and  one  introduced  by  this  author.  The  overall  goal  of  this 
research  was  to  analyze  algorithms  for  integer  factorization  based  on  the  use  of  continued 


fractions  and  quadratic  forms,  focusing  on  proving  mathematical  results  in  these  areas  and 
determining  average  runtime.  We  proposed  several  preliminary  sub-goals  to  this  end: 

1.  Analyze  the  conditions  for  which  Square  Forms  Factorization  provides  a  factorization. 

2.  Analyze  the  connection  between  continued  fractions  and  quadratic  forms  and  provide 
related  proofs. 

3.  Analyze  a  test  of  direction. 

4.  Produce  a  computer  implementation  of  these  algorithms. 

The  first  sub-goal  was  to  analyze  the  conditions  for  which  Square  Forms  Factorization 
provides  a  factorization.  This  information  is  important  to  cryptology,  as  the  security  of  some 
public  key  cryptosystems  is  dependent  on  the  difficulty  of  factorization.  The  conditions  for 
which  the  algorithms  must  provide  a  factorization  have  been  determined  but  a  sufficient  set 
of  conditions  for  which  the  algorithms  will  not  work  has  not  yet  been  determined.  How¬ 
ever,  the  runtimes  vary  greatly,  so  in  addition  to  the  original  goal,  this  research  analyzed 
which  numbers  Square  Forms  Factorization  may  factor  significantly  faster  or  slower  than 
others,  which  is  just  as  important  since  numbers  that  factor  extremely  slowly  are  almost  as 
secure  as  numbers  that  the  algorithm  doesn’t  factor  at  all.  It  is  possible  to  determine  suffi¬ 
cient  conditions  for  which  a  number  factors  quickly  but  necessary  conditions  remain  elusive. 
Specifically,  if  N  is  of  the  form  (a2m  +  b)(c2m  +  d),  it  is  possible  to  define  conditions  for 
special  cases,  but  these  conditions  have  not  been  generalized.  Proposition  1  of  §4.3  provides 
a  specific  example  of  this. 

The  second  sub-goal  was  to  analyze  the  connection  between  continued  fractions  and 
quadratic  forms,  specifically  intending  to  analyze  the  work  done  by  Shanks.  This  research 
has  taken  significantly  longer  than  planned  for  in  the  original  timeline.  Much  of  the  existing 
information  on  these  theories  was  scattered  and  disorganized,  and  many  of  the  papers  either 
lacked  proofs  or  contained  errors,  so  it  was  valuable  to  organize  this  information  into  a  form 
that  could  be  useful.  The  most  important  achievement  in  this  was  the  development  and 
proof  of  a  formula  relating  infrastructure  distance  with  composition  of  quadratic  forms. 
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The  third  sub-goal  was  to  analyze  a  test  of  direction  useful  to  a  new  algorithm  being  de¬ 
veloped.  The  research  in  this  area  has  demonstrated  that  a  test  of  direction  within  continued 
fractions  probably  cannot  be  used  effectively  for  factorization. 

The  fourth  of  the  original  sub-goals  was  to  produce  a  computer  implementation  of  these 
algorithms.  Shanks’  Square  Forms  Factorization  has  been  implemented  in  both  C  and 
Magma,  including  libraries  of  functions  concerning  the  operations  essential  to  quadratic 
forms  and  continued  fractions.  This  source  code  will  be  made  available  in  a  distributable 
form. 

Concerning  runtime,  Jason  Gower  [9]  has  recently  analyzed  the  runtime  of  SQUFOF 
for  his  Ph.D.  thesis  at  Purdue  University.  Included  in  this  was  an  analysis  of  multipliers, 
a  concept  that  will  be  used  in  these  implementations.  However,  it  has  been  conjectured 
by  Pomerance  that  as  SQUFOF  is  easy  to  parallelize,  a  parallel  implementation  may  be 
competitive,  so  a  new  goal  was  set  of  analyzing  the  parallel  implementation  of  SQUFOF. 

A  parallel  version  of  SQUFOF  has  been  developed  and  has  been  implemented  and  tested 
in  C.  This  code  is  included  in  Appendix  B.  This  algorithm  has  been  compared  with  an  older 
parallel  implementation  of  SQUFOF. 

Chapter  2  provides  some  background  to  the  problem  of  factorization,  basic  number  theory, 
and  some  of  the  tools  related  to  Square  Forms  Factorization.  Chapter  3  describes  how  Square 
Forms  Factorization  developed  from  the  existing  theory.  Chapter  4  describes  Square  Forms 
Factorization,  including  several  variations  of  the  algorithm,  and  describes  some  significance 
the  algorithm  may  have  for  cryptology. 
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Chapter  2 
Background 

2.1  The  Problem  and  its  Importance 

There  are  several  different  kinds  of  factorization;  this  research  will  focus  on  integer  factoriza¬ 
tion.  Consider  an  integer  N.  Factorization  is  the  process  of  finding  integers  p  and  q  greater 
than  1  such  that  N  =  pq.  Complete  factorization  would  require  repeating  this  process  for 
both  p  and  q  until  all  of  the  remaining  factors  are  prime.  However,  if  an  algorithm  can  be 
developed  to  quickly  factor  N  into  p  and  q,  the  same  algorithm  can  be  used  over  again  on  p 
and  q.  For  example,  it  is  easy  to  see  that  105  =  5-21  and  then  repeat  to  factor  21.  From 
here,  you  would  see  that  since  21  =  3  •  7,  105  =  5  •  3  •  7.  Although  in  this  simple  case,  the 
complete  factorization  is  easy  to  find,  this  task  becomes  much  harder  for  large  numbers. 

2.1.1  Number  Theoretic  Applications 

In  number  theory,  the  factors  of  a  number  dictate  many  of  the  characteristics  of  the  number. 
For  example,  Euler’s  0  function,  which  tells  how  many  numbers  less  than  N  are  relatively 
prime  to  N ,  can  be  directly  calculated  from  the  complete  factorization.  For  an  integer  that 
is  the  product  of  distinct  odd  primes,  (f>(N)  may  be  calculated  by  multiplying  each  of  the 
factors  minus  1. 

Example  1.  105  =  5  •  3  •  7.  Therefore  0(105)  =  (5  —  1)(3  —  1)(7  —  1)  =  48,  so  there  are  48 
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integers  less  than  105  that  are  relatively  prime  to  105. 

Also,  determining  whether  a  number  is  a  quadratic  residue  (i.e.  the  square  of  another 
number  modulo  N),  can  be  determined  directly  using  Gauss’s  Quadratic  Reciprocity  Law  if 
the  complete  factorization  of  N  is  known. 

Example  2.  192  =  361  =  46  +  3-105,  so  that  46  is  a  quadratic  residue  modulo  105  with  square 
root  19.  One  would  represent  this  as  192  =  46  (mod  105).  With  the  complete  factorization, 
we  can  use  Gauss’s  Quadratic  Reciprocity  Law  to  analyze  46  modulo  3,5,  and  7  to  determine 
whether  or  not  46  is  a  quadratic  residue  modulo  105  without  actually  having  to  find  its  square 
root  first  [8]. 

Therefore,  factoring  large  numbers  has  been  a  focus  of  research  for  a  variety  of  theoretical 
reasons. 

2.1.2  Cryptographic  Applications 

One  practical  application  of  factorization  is  public  key  cryptography.  The  idea  of  public  key 
cryptography  is  that  it  is  possible  to  keep  a  communication  secret  without  having  to  keep  the 
key  secret.  The  message  is  encrypted  by  one  key,  which  is  made  public,  and  is  decrypted  by 
another  key,  which  is  kept  secret.  In  the  most  prevalent  public  key  encryption  system,  RSA1 
[19],  the  cryptographer  chooses  a  number  that  is  the  product  of  two  large  prime2  numbers  p 
and  q:  N  =  pq.  Then  an  exponent  e  is  chosen  such  that  e  is  relatively  prime  to  (p—  l)(g  —  1). 
Although  there  are  several  variations,  in  the  normal  public  key  version,  N  and  e  are  made 
public.  The  user  of  the  RSA  system  then  privately  calculates  d,  the  inverse  of  e  modulo 
(p  —  l)(q  —  1).  Anyone  is  able  to  encrypt  something  to  him  by  raising  blocks  of  the  message 
to  the  power  e,  modulo  N:  c  =  me  (mod  N),  where  m  is  the  original  message  and  c  is  the 
encrypted  message.  Then,  the  recipient  is  able  to  decrypt  by  evaluating  m  =  cd  (mod  N) 

[19]- 

1RSA  is  named  after  Rivest,  Shamir,  and  Adleman.  It  was  earlier  developed  by  Clifford  Cooks  of  GCHQ, 
but  this  was  only  recently  declassified  [5]. 

2Usually,  these  are  just  numbers  that  pass  several  primality  tests  and  thus  have  a  high  probability  of 
being  prime,  called  pseudo-primes.  Very  rarely,  one  of  them  will  not  be,  but  this  is  rare  enough  to  not  cause 
significant  problems. 
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Example  3. 


p  =  4327,  q  =  1009,  N  =  pq  =  4365943 
e  =  2005 


and  from  this  the  cryptographer  would  use  p,  q  and  e  to  calculate  d  =  865597.  He  would 
make  N  and  e  public.  Suppose  someone  wants  to  send  him  the  message  809. 

The  sender  would  evaluate  80  92005  (mod  N)  as  1591327.  This  would  be  sent  over  the 
internet. 

The  receiving  computer  would  then  evaluate  159  132786559'  (mod  N)  as  809.  Since  this  is 
the  only  computer  with  access  to  d,  anyone  else  intercepting  the  message  would  have  great 
difficulty  in  determining  the  message 

Another  application  of  public  key  cryptography  is  for  signatures.  If  the  sender  uses 
his  private  key  to  encrypt  a  message,  the  receiver  may  use  the  public  key  to  decrypt  the 
message  and  be  certain  of  who  the  message  originated  from  and  that  no  one  along  the  way 
has  modified  it.  Often  this  is  combined  with  a  second  layer  of  encryption  to  prevent  anyone 
else  from  also  being  able  to  use  the  public  key  to  read  the  message. 

As  long  as  someone  intercepting  a  message  is  unable  to  factor  N,  it  is  usually  impossible 
to  obtain  d,  so  that  the  message  cannot  be  broken.  The  security  of  RSA  and  its  variations 
depends  highly  on  whether  or  not  N  can  be  factored  [24],  Although  extremely  fast  factor¬ 
ization  would  be  a  threat  to  these  systems,  the  advances  in  number  theory  produced  by 
faster  factorization  would  likely  provide  a  number  of  alternative  secure  systems.  Also,  in 
addition  to  the  potential  for  alternative  secure  systems,  there  is  the  strong  possibility  that 
fast  factorization  algorithms  will  work  better  for  some  numbers  than  others.  If  there  are 
classes  of  numbers  that  a  faster  factorization  algorithm  does  not  work  on,  this  would  en¬ 
able  designers  of  the  algorithm  to  increase  their  security  by  relying  more  on  these  numbers. 
Regardless  of  whether  or  not  the  algorithm  works  for  all  numbers  or  provides  alternative 
systems,  for  security  purposes  it  is  necessary  to  understand  the  strengths  and  weaknesses  of 
the  algorithm. 
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2.1.3  Runtime  Analysis 

Up  to  this  point  we  have  referred  to  the  speed  of  factorization  in  general  terms,  but  there  are 
several  different  ways  to  classify  the  speed  of  an  algorithm.  Let  N  be  the  number  to  factor. 
Let  n  be  the  number  of  bits  in  N,  n  =  log2  N.  An  algorithm’s  run  time  is  called  exponential 
if  it  increases  exponentially  with  the  size  of  the  input,  in  this  case  n.  Linear  refers  to  an 
algorithm  where  the  time  increases  proportionally  to  the  number  of  bits3.  Polynomial  refers 
to  an  algorithm  for  which  the  time  required  is  some  polynomial  function  of  n.  Thus,  linear 
time  is  a  special  case  of  polynomial  time.  There  are  some  algorithms  that  fall  in  between 
polynomial  and  exponential  time  and  are  referred  to  as  sub-exponential. 

Note  that  different  algorithms  may  be  faster  for  different  sizes  of  numbers  or  even  for 
different  systems.  This  is  why  the  runtime  is  analyzed  in  terms  of  growth.  For  small  values 
of  n.  linear  and  exponential  may  be  fairly  close  or  linear  may  even  be  faster,  but  the  runtime 
of  an  exponential  algorithm  will  grow  faster  with  the  size  of  n,  so  that  for  sufficiently  large 
n,  an  exponential  algorithm  will  be  slower  than  a  sub-exponential,  which  will  be  faster  than 
a  polynomial,  which  will  be  faster  than  a  linear. 

The  runtimes  for  sub-exponential  algorithms  are  described  by  complicated  formulas.  De¬ 
fine 

L(a,c )  =  exp(cn"(logn)1'")) 

If  a  =  1,  this  is  an  exponential  algorithm.  If  a  =  0,  this  is  a  polynomial  algorithm.  For 
values  of  a  in  between  0  and  1,  the  algorithms  are  sub-exponential. 

In  terms  of  asymptotic  average  runtime,  the  best  general  purpose  factorization  algorithm 
is  the  general  number  field  sieve ,  with  a  runtime4  of  L(  1/3, 7^73  +  o(l))  [6].  Two  other 
algorithms  are  currently  in  common  use,  the  elliptic  curve  method ,  which  has  a  runtime  of 
L(l/2,1  +  o(l))  [13]  and  the  multiple  polynomial  quadratic  sieve ,  which  has  a  runtime  of 

3Since  n  =  log2  N,  so  that  such  a  runtime  is  logarithmic  in  N,  this  is  often  referred  to  as  logarithmic, 
resulting  in  a  certain  amount  of  confusion.  Note  that  both  ways  of  expressing  runtime  will  be  used  in  this 
paper. 

4The  expression  o(l)  is  little  o  notation.  It  converges  to  0  as  n  goes  to  infinity.  An  understanding  of  little 
o  notation  is  not  important  for  understanding  this  report. 
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L(  1/2, 1)  [16].  Although  each  of  these  algorithms  is  slower  on  average,  they  are  each  faster 
for  some  types  of  integers,  so  that  in  combination  they  tend  to  be  faster  than  the  number 
field  sieve. 

Another  less  commonly  considered  characteristic  of  a  factorization  algorithm  is  the  effi¬ 
ciency  with  which  it  can  use  multiple  processors  working  simultaneously.  The  current  trend 
in  ultra-fast  computing  is  to  use  a  large  number  of  inexpensive  processors  instead  of  a  single 
expensive  processor.  Given  this  trend,  it  is  important  that  a  good  algorithm  be  able  to  use 
a  multiple  processor  system  efficiently.  Dividing  a  single  task  between  multiple  processors  is 
called  parallelization.  With  perfect  efficiency,  factoring  a  number  with  10  processors  should 
take  on  average  one  tenth  the  time  required  for  factorization  on  a  single  processor,  but  no 
algorithm  can  ever  quite  attain  perfect  efficiency. 

2.2  Basic  Number  Theory 

Mathematics  is  the  Queen  of  the  sciences,  and  arithmetic  the  Queen  of  mathe¬ 
matics. 


C.  F.  Gauss  [25] 

This  section  provides  a  quick  overview  of  basic  number  theory  concepts  and  notation 
that  will  be  used  throughout  this  report.  Those  already  familiar  with  number  theory  are 
welcome  to  skip  this  section,  as  none  of  this  is  original  and  may  be  found  in  any  standard 
number  theory  text  [10]. 

Number  theory  analyzes  sets  of  numbers.  If  an  element  e  is  in  a  set  S,  this  is  written 
e  G  S.  If  an  element  e  is  not  in  a  set  S,  this  is  written  e  ^  S'.  If  all  the  elements  in  set  Si  are 
in  set  S2,  then  Si  is  contained  in  set  S2  and  this  is  written  Si  C  S2.  If  Si  C  S2  and  S2  C  Si, 
then  Si  =  S2. 

The  set  of  all  real  numbers  is  denoted  M.  Although  this  set  will  be  important,  it  is 
primarily  several  special  subsets  of  the  real  numbers  that  this  research  will  focus  on: 
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The  first  subset  of  the  real  numbers  is  the  set  of  integers,  represented  as  Z. 

Z  =  {...-3,-2,-l,0,l,2,3...}. 

The  natural  numbers  are  the  integers  >  1.  The  operations  on  the  integers  are  the  standard 
addition  and  multiplication.  Addition  has  the  usual  four  basic  properties:  associativity, 
commutativity,  identity,  and  inverses. 

Multiplication  has  three  basic  properties:  associativity,  commutativity,  and  identity. 

One  other  property,  distributivity,  involves  both  addition  and  multiplication: 
a  ■  (b  +  c)  =  ab  +  ac. 

There  are  several  other  relations  that  are  important  in  Z.  A  natural  number  a  divides 
an  integer  b  if  there  exists  an  integer  k  such  that  b  =  ka  and  this  is  written  a  \  b.  If  a  does 
not  divide  b,  this  is  written  a  \  b.  A  natural  number  >  1  is  prime  if  it  may  only  be  divided 
by  itself  and  1.  There  are  several  basic  properties  of  this  relation: 


If  a 

b  and  b 

a,  then  a  - 

=  b. 

If  a 

b  and  a 

c ,  then  a 

( b  +  c 

If  a 

b,  then  a 

be. 

If  a  is  prime  and  a  \  be,  then  a  \  b  or  a  \  c. 

If  for  some  integer  n  and  a  prime  p,  pn  \  b,  but  pn+1  \  b,  this  is  written  pn  ||  b. 
c  is  the  greatest  common  divisor  of  a  and  b  if  c  >  0,  c  |  a,  c  \  b,  and  for  any  integer  d  such 
that  d  |  a  and  d  \  b,  d  \  c.  This  is  written  c  =  gcd(a,  b)  or  sometimes  (although  not  in  this 
report)  as  merely  (a,  b). 

Two  integers  a  and  b  are  relatively  prime  if  gcd(a,  b)  =  1,  that  is,  if  a  and  b  have  no  factor 
in  common. 

If  p,  m,  n  are  integers  with  p  prime  such  that  pm  ||  a  and  pn  ||  b,  then  pmm(m-n)  || 

ged  (a,b). 

If  c  =  gcd(a,  b),  there  exist  integers  x,  y  G  Z  such  that  ax  +  by  =  c.  gcd(a,  b),  along  with 
these  two  integers,  may  be  efficiently  calculated  using  Euclid’s  extended  algorithm. 

c  is  the  least  common  multiple  of  a  and  b  if  c  >  0,  a  \  c,  b  |  c,  and  for  any  integer  d  such 
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that  a  |  d  and  b  \  d,  c  \  d.  This  is  written  c  =  lcm(a,  b)  or  sometimes  merely  as  {a,b}.  If 
p,m,n  are  integers  with  p  prime  such  that  pm  ||  a  and  pn  ||  b,  then  pmax(m>n)  ||  lcm(a,5). 

The  gcd  and  1cm  are  related  by  the  formula  gcd(a,  b)  ■  lcm(a,  b)  =  ab. 

If  c  j  ( a  —  b ),  then  we  say  that  a  is  congruent  to  b  modulo  c,  written  a  =  b  (mod  c).  This 
relation  has  several  basic  properties: 

If  a  =  b  (mod  n),  then  b  =  a  (mod  n). 

If  a  =  b  (mod  n)  and  b  =  c  (mod  n),  then  a  =  c  (mod  n). 

If  a  =  b  (mod  n)  and  c  =  d  (mod  n),  then  a  +  b  =  c  +  d  (mod  n)  and  ac  =  bd 

(mod  n) . 

By  using  Euclid’s  extended  algorithm  to  compute  greatest  common  divisors,  for  a  rela¬ 
tively  prime  to  b  it  is  possible  to  calculate  ar1  (mod  b).  Thus,  division  modulo  b  is  defined 
as  reducing  to  least  terms  and  then  multiplying  by  the  inverse. 

One  important  type  of  number  is  a  quadratic  residue.  The  integer  a  is  a  quadratic  residue 
of  integer  base  b  if  there  exists  some  integer  x  such  that  x2  =  a  (mod  b).  For  example,  2  is 
a  quadratic  residue  of  7  because  32  =  9  =  2  (mod  7).  The  symbol  used  is  (|).  (|)  =  1  if 
x2  =  a  (mod  b)  has  a  solution,  —1  if  it  does  not,  and  0  if  a  and  b  are  not  relatively  prime.  If 
a  is  a  quadratic  residue  of  b  and  a  is  a  quadratic  residue  of  c  and  b  and  c  are  relatively  prime, 
then  a  is  a  quadratic  residue  of  be.  Gauss’s  Quadratic  Reciprocity  Law  provides  a  fast  way 
of  determining  quadratic  residues  modulo  a  base  whose  prime  factorization  is  known: 

Theorem  2.2.1.  For  p  and  q  distinct  primes: 

( ^  (5)  =  (_!)(P— 1)(9— 1)/4 
\qj  V 

(2)  = 

(y)  =  (-i)(r-1,/2 

There  are  many  other  well-known  properties  concerning  congruences.  For  example,  “Fer¬ 
mat’s  little  theorem”:  for  p  prime  and  a  an  integer  if  p  \  a,  then  ap_1  =  1  (mod  p).  This 
is  often  used  to  test  whether  not  an  integer  N  is  prime.  (If  for  some  a  <  N,  aw_1  ^  1 
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(mod  N),  then  N  is  not  prime.  However,  there  are  infinitely  many  composite  numbers  that 
will  appear  prime  for  any  a  such  that  gcd (a,N)  =  1,  so  the  converse  is  not  necessarily  true 

[l]-) 

In  Z,  inverses  are  not  defined  for  multiplication.  However,  this  question  brings  us  to  the 
set  of  rational  numbers,  represented  as  Q. 

Q  =  {a/b  :  a,  b  G  Z,  b  ^  0} 

All  of  the  properties  for  addition  and  multiplication  in  Z  still  apply,  in  addition  to  one 
more: 

If  a  G  Q  and  a  ^  0,  then  there  exists  a  unique  integer  b  such  that  a  ■  b  =  1.  This 
is  written  as  aG1  or  1/a. 

Let  a  denote  a  real  or  complex  number  not  belonging  to  Q.  Denote  the  set  of  all  finite 
linear  combinations  with  integer  coefficients  of  all  the  non-negative  powers  of  a  by, 


Z[a]  —  T  CL\Ot  T  02*3:^  T  030:^  T  ...  :  a*  G  Zj-. 

For  example,  if  a  =  \/2,  then  Z[\/2]  contains  5  +  3^  and  4  —  2-\/2,  but  does  not  contain 
Alternately,  the  new  set  could  be  the  set  of  all  linear  combinations  of  all  the  powers  of 
a.  That  is, 


Q(o:)  —  {...a_2Q:  ^  T  (X—\Ot  '  T  clq  T  cl\Ql  T  0,2^  T  •••  :  a.j  G  Q} 

For  example,  Q(v^2)  contains  both  3  +  \/2  and  4  +  ^,  but  does  not  contain  y^. 

For  an  element  (  in  either  of  the  sets  Z[VfV]  and  Q (x/fV),  (  refers  to  the  conjugate  of  </, 
which  is  found  by  changing  the  sign  of  the  algebraic  part.  For  example,  if  (  =  1  +  v^3,  then 

c  =  1 -  V3. 

The  norm  of  (  is  J\f( Q  = 
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2.3  Quadratic  Field  Number  Theory 

This  section  introduces  four  areas  of  mathematics  that  this  research  has  focused  on:  con¬ 
tinued  fractions,  quadratic  forms,  ideals,  and  lattices.  These  areas  are  closely  related  and 
therefore  have  the  same  property  of  being  periodic.  See  the  appendices  for  a  more  detailed 
analysis  of  each  of  these  topics,  along  with  related  proofs  and  an  analysis  of  their  connections. 

1.  A  continued  fraction  is  a  tool  originally  used  for  rational  approximation.  Given  a 
number  a,  the  goal  is  to  represent  a  in  the  form 


OL  —  6q  + 


1 


6i  + 


i 

62  +  ... 


for  integers  6*.  This  is  often  abbreviated  as  [60,  b\,  b-2, ...].  The  expressions  found  by 
truncating  this,  b0,  [60,  ffi],  [&o?  &i,  £>2],  •••  are  called  the  convergents.  These  simplify  down 
to  rational  numbers  An/Bn. 


The  sequence  of  bf  s  is  found  by  the  recursive  formulas 


x0  =  a,  b0=  |x0J  (2.1) 

Vi  >  1  Xi  =  - — — ,  bi  =  [xi\  (2.2) 

Xi-i  -  bi- 1 

Note  that  this  formula  is  derived  by  solving  the  equation  1  =  bi- 1  +  A  for  ay. 
refers  to  the  floor  of  x,  the  greatest  integer  less  than  or  equal  to  x. 

If  Xfl  =  \^N  for  N  a  non-square  integer  (that  is,  not  the  square  of  another  integer),  then 
this  sequence  may  be  calculated  efficiently.  In  particular,  each  ay  reduces  to  the  form 
A) (+p  with  P  and  Q  integers.  The  integer  P  in  the  numerator  is  called  the  residue. 
The  denominator  Q  is  called  a  pseudo-square.  There  are  formulas  to  calculate  these 
recursively,  after  the  first  several  steps,  without  doing  arithmetic  on  any  integer  larger 
than  yfN : 
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k 


[Vn\  +  Pj_ i 

Qi 


Pi  —  biQi  —  Pi-i  Qi+i  —  Qi- 1  +  bi(Pi_  i  —  Pi) 


(2.3) 


The  rational  approximation  produced  may  also  be  evaluated  recursively,  so  that  the 
integers  and  Bi  such  that  Ai/Bi  =  [60,  &i,  -  --&*]  may  be  calculated  efficiently.  For 
factorization,  the  relation 


Al^i-lfQi  (mod  N) 


(2.4) 


has  been  important  for  several  algorithms.  This  equation  also  explains  the  name 
‘pseudo-square’. 

Take  the  following  example: 

Example  4. 


Xl  = 

x2  = 

X3  = 


x0  = 

1 T— ^ 

II 

b0  =  6 

1 

i 

VS+6 

1 

1 

o 

a 

V41-6 

5 

i 

5 

v^I+4 

xi-b  i 

FiI-4 

5 

1 

5 

vAT+6 

x2-b2 

FiI-6 

1  ’ 

>  evident  that 

X4  =  X\ 

VAL  =  6  +  d- 


,  h  =  2  y/41  =  6  + 

,  b2  =  2  VU  =  6  +  - 


Xl 


2+  — 

x2 


2+- 


2+- 


occurs  with  the  continued  fraction  expansion  for  \/]V  for  N  non-square.  The  length  of 
this  period  is  closely  related  to  several  deep  number  theory  ideas  involving  the  class 
number  and  the  regulator. 


A  more  thorough  explanation  of  continued  fractions,  including  detailed  proofs  is  given 
in  Appendix  A.l. 

2.  A  binary  quadratic  form  is  a  polynomial  of  the  form  F(x,y )  =  ax2  +  bxy  +  cy2,  for 
a,  b,  c  G  Z.  (Often  this  is  abbreviated5  as  (a,  b,c )). 


5 As  b  is  usually  even,  many  writings  on  quadratic  forms  represent  this  as  the  triplet  (a,  6/2,  c),  so  that 
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Two  quadratic  forms  are  equivalent  if  they  have  the  same  range,  the  set  of  possible 
values  for  F(x,y)  for  integer  values  of  x  and  y,  written  F\  ~  F2.  Appendix  A. 2  gives 
a  more  precise  description  of  equivalence.  The  surprising  fact  is  that  two  quadratic 
forms  are  equivalent  if  and  only  if  they  correspond  to  terms  from  the  same  continued 
fraction  expansion  (Theorem  A. 2. 4  in  Appendix  A. 2),  so  that  the  cycle  formed  by  the 
XiS  used  to  compute  continued  fractions  corresponds  to  a  cycle  of  equivalent  quadratic 
forms. 

The  primary  operation  on  quadratic  forms  is  composition.  This  operation  is  defined  in 
detail  in  Appendix  A. 2,  Proposition  2.  The  composition  of  two  binary  quadratic  forms 
is  another  binary  quadratic  form,  written  with  the  symbol  (*).  Note  that  although 
the  original  definition  of  this  operation  is  related  to  multiplication,  this  is  an  operation 
that  is  quite  distinct  from  multiplication,  as  the  product  of  two  binary  quadratic  forms 
is  no  longer  a  binary  quadratic  form. 

Equivalence  and  composition  are  closely  related.  Specifically,  if  Fi  ~  F2,  then  F2  *G  ~ 
F2  *  G.  The  cycles  of  equivalent  forms  are  the  elements  of  the  class  group,  which  has 
been  studied  extensively. 

A  more  thorough  explanation  of  quadratic  forms  is  given  in  Appendix  A. 2. 

3.  A  lattice  is  the  set  of  all  finite  linear  combinations  (with  integer  coefficients)  of  a 
generating  set  that  contains  a  vector  space  basis  for  the  rational  vector  space  the  lattice 
is  contained  in.  For  this  research,  the  ambient  vector  space  is  Q(-\/A0  x  Q(VN).  If 
OL\ 5  0^2 1  •  •  •  5  eQ  (Vn)  x  Q(v^V)  are  vectors,  then 

(  k 

,  OL2,  ...,  Ctfc]  \  ^  ^  'TliOLi  .  Hi  £  Zi 

l  i=  1 

is  a  lattice. 

The  notation  is  often  abused  slightly,  so  that  the  vector  {(,  ()  is  represented  by  just  (. 


— 14x2  + 10 xy+5y2  would  be  represented  (—14,  5,  5).  This  report  does  not  take  out  this  factor  of  2.  Although 
this  will  be  clear  throughout  this  report,  please  observe  the  potential  for  confusion  on  this  “annoying  little 
2  [2]”  when  comparing  different  sources  on  quadratic  forms. 
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There  are  several  properties  of  lattices  that  are  important  for  this  research. 
Definition  2.3.1.  For  a  vector  v  =  (i’i,v2),  the  normed  body  of  v,  7 Z(v)  is  the  set 

TZ{y)  =  {(a?i, a?2)  :  x\,x2  e  M,  \x±\  <  |ui|,  \x2\  <  |u2|} 

Abusing  notation  again,  denote  7 Z(£)  =  7£((£,  £)). 

A  number  £  (or  actually  the  corresponding  vector)  is  a  minimum  of  £  if  7£(£)fl£  =  {0}, 
where  0  is  the  vector  (0,0).  Figure  2.1  provides  an  example  of  this. 

A  lattice  £  is  reduced  if  1  e  £  and  1  is  a  minimum. 


• 

u 

w 

y 

Figure  2.1:  A  Lattice  with  Minima 


The  image  above  demonstrates  minima  of  a  lattice.  The  lattice  is  the  set  of  points. 
Each  of  the  points  u,  v,  and  w  are  minima  of  this  lattice,  since  there  are  no  points  of 
the  lattice  inside  their  respective  normed  bodies  (the  rectangles)  except  the  origin. 

Another  important  characteristic  of  minima  is  adjacency.  If  (xi,yi)  and  (x2,  y2)  are 
minima  with  |xi|  >  \x2\  and  \y\  <  y2 \ ,  these  two  minima  are  adjacent  if  there  does 
not  exist  another  minima  (x3,y3)  such  that  \x2\  <  |x3|  <  |o?i |  and  \yi\  <  \y3\  <  \y2\.  In 
Figure  2.1,  u  and  v  are  adjacent  minima  and  v  and  w  are  adjacent  minima. 

The  important  connection  between  lattices  and  continued  fractions  is  a  method  devel¬ 
oped  by  Voronoi  for  finding  a  chain  of  adjacent  minima  for  a  lattice  (and  a  chain  of 
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other  reduced  lattices).  This  process  matches  the  continued  fraction  algorithm  exactly, 
a  property  that  implies  that  the  distance  derived  from  lattices  also  applies  to  continued 
fractions  and  quadratic  forms. 

A  more  thorough  description  of  lattices  is  given  in  Appendix  A. 4. 

4.  An  ideal  is  a  common  algebra  concept  that  is  used  in  a  variety  of  contexts.  For  a  ring 
S'  ,  a  set  /  is  an  ideal  of  S  if 

I  CS. 

If  a,  b  E  I,  then  a  +  b  E  I  and  a  —  b  E  I. 

If  a  E  I  and  s  E  S,  then  a  ■  s  E  I. 

A  classic  example  is  the  set  of  all  even  integers,  2Z,  within  the  set  Z  of  all  integers. 
The  sum  of  two  even  integers  is  still  an  even  integer  and  the  product  of  an  even  integer 
with  any  integer  is  an  even  integer. 

This  analysis  focuses  on  ideals  of  Z  [VA/V],  for  N  a  non-square  positive  integer.  Describ¬ 
ing  these  ideals  will  again  require  the  notation  for  the  lattice  generated  by  a  set.  With 
this  notation  the  ideals  of  Z [\/N]  are  of  the  form  [Q,  sy/N +P],  for  some  integers  Q ,  s,  P 
and  a  set  I  =  [Q,s\/N  +  P]  is  an  ideal  of  Z [y/N]  if  and  only  if  sQ  |  A f(s\/N  +  P), 
s  j  Q ,  and  s  \  P  (Appendix  A. 3,  Lemma  A. 3. 2).  An  ideal  is  primitive  if  s  —  1. 

The  multiplication  of  ideals  is  important  to  many  fields.  For  the  ideals  /  =  [Q,  \/N +P] 
and  J  =  [Q',y/N  +  P'], 

I  ■  J  —  [QQ\  qVn  +  QP\  Q'Vn  +  Q'P,(Vn  +  P)(Vn  +  P)} 

This  may  then  be  simplified  to  reduce  it  to  an  ideal  of  the  form  s[q,  y/N  +p\,  where  s, 
q,  and  p  may  be  calculated  efficiently. 
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The  quadratic  form  (of  discriminant  A  =  0  (mod  4))  F(x,y)  =  Ax 2  +  Bxy  +  Cy 2 
corresponds  to  the  ideal 


A, 


B 

~2 


and  similarly  the  ideal  [Q,  +  P]  corresponds  to  the  quadratic  form 


F(x,  y)  =  Qx 2  +  2 Pxy  + 


P2  —  N 

Q 


Note  that  A  =  AN. 


The  value  of  ideals  is  that  the  equations  defining  multiplication  of  ideals  correspond 
exactly  to  the  equation  defining  composition  of  quadratic  forms,  providing  a  connection 
between  quadratic  forms  and  lattices.  A  more  thorough  description  of  ideals  is  given 
in  Appendix  A. 3. 


The  appendices  also  describe  the  one  to  one  correspondence  between  the  elements  of  these 
different  sets.  For  this  paper,  it  is  sufficient  to  understand  that  these  maps  do  exist,  that 
is,  each  term  in  the  continued  fraction  expansion  is  paired  with  exactly  one  quadratic  form, 
exactly  one  lattice,  and  exactly  one  ideal  and  given  any  one,  any  of  these  four  expressions 
may  be  calculated  immediately.  The  indices  derived  from  continued  fractions  are  typically 
also  used  for  the  related  quadratic  forms,  ideals,  and  lattices. 

As  all  of  these  tools  are  closely  related,  they  all  have  the  property  of  being  periodic,  as 
demonstrated  by  Example  4.  The  reason  for  this  is  that  the  bounds  on  the  coefficients  of  a 
reduced  quadratic  form  (Definition  A. 2.1)  imply  that  there  can  only  be  a  finite  number  of 
them.  Since  the  recursive  formulas  (2.3)  provide  a  way  of  getting  from  one  reduced  quadratic 
form  to  the  next,  this  process  must  eventually  repeat  and  thus  be  periodic.  The  length  of 
this  period  is  equal  to  the  regulator  of  the  class  group  and  is  related  to  many  other  number 
theory  concepts  ([12],  [20]). 
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If  the  continued  fraction  expansion  is  started  with  something  other  than  xq  =  V7V,  the 
result  may  be  a  sequence  that  is  completely  disjoint  from  the  principal  cycle  (the  sequence 
that  begins  with  xq  =  \/N).  Most  of  these  cycles  are  symmetric6.  Lemma  A.  1.6  and  Theorem 
A. 1.7  provide  a  thorough  definition  and  analysis  of  these  symmetries,  but  the  concept  can 
be  readily  understood  from  several  examples: 

Example  5.  In  Example  4,  the  period  of  the  denominators  was  1,5,5, 1, ...,  so  that  1  and  5 
are  both  symmetry  points.  If  Xo  =  a/115,  then  the  period  of  the  denominators  is 

1,15,6,11,9,10,9,11,6,15,1,... 

so  that  the  symmetry  points  are  1  and  10. 

If  instead  Xq  =  v/^+9,  then  the  period  of  the  denominators  is 

2,17,3,5,3,17,2,... 

This  is  a  different  cycle  for  the  class  group  for  N  =  115. 

These  are  the  cycles  that  are  important  for  factorization. 

The  square  of  a  quadratic  form  in  the  principal  cycle  is  also  in  the  principal  cycle  and 
Theorem  A. 2. 8  shows  that  the  square  of  a  form  in  any  ambiguous  cycle  is  also  in  the  principle 
cycle.  This  property  enables  a  variety  of  methods  of  finding  points  in  ambiguous  cycles  other 
than  the  principal  cycle. 


6 The  non-symmetric  cycles  are  not  analyzed  in  this  paper. 
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Chapter  3 

History:  From  Morrison  and  Brillhart 
to  Present 

3.1  Before  Shanks 

In  1931,  D.  H.  Lehmer  and  R.  E.  Powers  [11]  introduced  one  simple  and  fairly  intuitive 
algorithm  for  using  continued  fractions  for  factorization.  Unfortunately,  as  it  had  the  frus¬ 
trating  tendency  to  fail  a  few  times  before  succeeding,  it  was  dropped.  However,  as  better 
computing  capabilities  became  available,  John  Brillhart  began  to  consider  that  although 
this  algorithm  was  cumbersome  for  smaller  numbers,  it  might  have  an  advantage  for  larger 
numbers.  In  1970,  Morrison  and  Brillhart  implemented  CFRAC,  as  it  became  known,  and 
tried  to  factor  a  39  digit  number  [15].  The  entire  algorithm  is  a  bit  more  complicated,  but 
here  is  a  description  sufficient  for  our  purposes. 

If  the  equation 


x2  =  y2  (mod  N) 

can  be  solved  such  that 


(3.1) 


x  ^  ±y  (mod  N ), 


(3.2) 
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then  gcd (x  —  y,  N)  provides  a  nontrivial  factor  of  N.  Choosing  a  value  for  x  and  then  looking 
for  a  value  for  y  that  satisfies  equations  (3.1)  and  (3.2)  is  not  computationally  effective  for 
large  N.  Fermat,  the  first  to  employ  this  concept,  tested  values  of  x  greater  than  a /N  to 
find  a  value  such  that  x2  (mod  N)  was  already  a  perfect  square.  However,  these  numbers 
get  large  quickly.  Continued  fractions  provide  another  approach  to  achieving  this  and  since 
0  <  Qi  <  2 '/N,  the  chances  of  finding  a  perfect  square  are  greatly  improved.  From  (2.4), 
Af_1  =  (— 1  YQi  (mod  N).  If  for  some  i,  Q21  is  a  perfect  square  then  this  provides  a  solution 
to  (3.1)  and  it  only  remains  to  check  whether  or  not  gcd(x  +  y,  N )  or  gcd(x  —  y ,  N )  provide 
a  nontrivial  factor  of  N.  However,  there  are  not  very  many  perfect  squares  in  the  continued 
fraction  expansion  so  Morrison  and  Brillhart  [15]  used  products  to  obtain  squares.  For 
example,  for  N  =  1333,  the  continued  fraction  expansion  provides  732  =  —3  (mod  N )  and 
17892  =  —12  (mod  N).  From  this,  (73-1789)2  =  (— 3)(— 12)  =  62  (mod  N),  quickly  yielding 
gcd(73  •  1789  -  6, 1333)  =  43,  so  that  1333  =  31  •  43. 

There  are  two  problems  with  this  algorithm.  First,  it  requires  the  calculation  of  the 
Hj’s,  which  are  of  the  same  size  as  N,  after  reduction  modulo  N,  while  the  rest  of  the 
algorithm  only  requires  arithmetic  on  numbers  of  size  a /N.  Second,  after  going  through  a 
nontrivial  amount  of  computation  to  find  a  relation  that  solves  (3.1),  not  all  of  these  result 
in  a  factorization.  I  provide  one  example: 

Example  6.  In  the  continued  fraction  for  ^1333,  Qe  =  9  =  32  and  A5  =  10661,  so  that 
106612  =  32  (mod  1333).  Unfortunately  10661  =  —3  (mod  1333),  so  this  square  does  not 
result  in  a  nontrivial  factor  of  1333. 


3.2  Shanks’  Work 

In  1975,  Daniel  Shanks  developed  several  very  interesting  algorithms  for  factorization  ([22], [23]) 
based  on  an  understanding  of  quadratic  forms  and  the  class  group  infrastructure.  First  he  de¬ 
veloped  an  improvement  to  the  Morrison-Brillhart  algorithm.  Roughly  speaking,  rather  then 
saving  the  Hj’s,  he  was  able  to  use  composition  of  quadratic  forms  to  combine  numbers  to 
produce  squares  and  then  use  the  “infrastructure”  to  use  those  squares  to  find  a  factorization. 
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In  addition,  he  developed  from  the  concept  of  infrastructure  a  system  of  predicting  whether 
or  not  any  given  square  would  provide  a  nontrivial  factor.  Unfortunately,  this  didn’t  save 
very  much  time  and  was  a  much  more  complicated  algorithm  than  the  Morrison-Brillhart 
algorithm. 

From  here  the  development  of  the  algorithm  was  prompted  by  the  number  260  +  230  — 1.  It 
failed  a  Fermat  primality  test1,  but  when  Morrison  and  Brillhart  tried  to  factor  it,  it  failed  114 
times.  Therefore,  they  stopped,  multiplied  it  by  some  small  constant  and  tried  again.  This 
time  it  worked  on  the  first  try,  but  they  wanted  to  know  why  it  had  failed  so  many  times. 
So  they  asked  Shanks  to  analyze  it.  Unfortunately  (or  fortunately  in  hindsight),  Shanks 
only  had  an  HP-65  available  and  he  couldn’t  fit  his  entire  algorithm  into  it.  Therefore,  he 
discarded  all  the  work  of  combining  numbers  to  form  squares  and  just  cycled  through  until 
he  found  one  already  there.  The  code  for  this  was  much  shorter,  and  as  it  turned  out  the 
algorithm,  which  became  known  as  Square  Forms  Factorization  or  SQUFOF,  was  actually 
significantly  faster.  Here’s  how  it  works: 

The  continued  fraction  cycles  contain  symmetry  points  that  potentially  provide  a  factor¬ 
ization  for  N  (Lemmas  A.1.6-A.1.7).  The  goal  of  SQUFOF  and  any  related  algorithms  is 
to  search  for  these  symmetry  points.  SQUFOF,  described  in  Figure  3.1,  provides  a  fairly 
simple  way  of  achieving  this.  The  basic  version  searches  the  sequence  in  order  until  it  finds  a 
perfect  square  on  an  even  index2.  Then  it  takes  the  square  root  of  the  denominator  and  the 
conjugate  of  the  numerator.  According  to  Theorem  A. 2. 8,  this  new  term  still  corresponds 
to  an  ambiguous  cycle,  that  is,  it  will  still  have  two  symmetry  points.  By  Theorem  A. 5. 2, 
one  of  these  symmetries  point  is  behind  this  term  (in  front  of  it  after  taking  the  conjugate) 
a  distance  that  is  half  of  the  distance  from  the  start  to  the  square  found.  Therefore,  the  al¬ 
gorithm  continues  cycling  from  here  until  it  finds  the  symmetry  point,  which  is  indicated  by 
a  repeated  numerator.  At  this  point  P*  must  then  have  a  factor  in  common  with  Qi  because 

1The  Fermat  primality  test  is  based  on  Fermat’s  classical  result  that  for  p  prime,  ap  =  a  (mod  p)  [10]. 
Equivalently,  if  aN  ^  a  (mod  N)  for  some  integer  a,  then  N  isn’t  prime. 

2 At  worst  case,  the  1  at  the  end  of  the  cycle  will  provide  this.  The  analysis  by  Shanks  [23]  and  later  by 
Gowers  [9]  predict  that  for  N  with  k  +  1  distinct  factors,  there  should  be  on  average  ln(8)  steps 

until  the  first  square  that  provides  a  non-trivial  factor. 


of  the  formula  Pi  =  biQi  —  Pt-\  (see  Case  1  of  Theorem  A. 1.7).  Since  N  =  Pf  +  QiQi+ 1, 
this  common  factor  is  a  factor  of  2N  and  may  be  found  by  taking  the  greatest  common 
divisor  of  P%  and  Qi.  For  an  analysis  of  how  Shanks  predicted  when  the  factor  found  would 
be  non-trivial,  see  §A.6.  With  the  inclusion  of  an  idea  Shanks  called  Fast  Return,  explained 
shortly,  these  predictions  will  be  unnecessary. 


Given  N  to  be  factored: 
x0  <-  VN  b0  <-  |_£0J 
while  Qi  ^  perfect  square 
Apply  equation(2.2)  twice 
Reduce  xt  to  the  form  '^0  Pi 

/  v (N+Pj 

VQl 

while  Pi  ^  Pi_i 

Apply  equation  (2.2)  and  reduce 

gcd (Qi,Pi),  the  greatest  common 
divisor,  is  a  nontrivial  factor  of  N. 

Figure  3.1:  SQUFOF  without  Composition 


Example  7.  Let  N  be  1353: 

x0  =  a/1353  b0  =  36 

1  _  a/I353  +  36  _  VT353  -  21 

Xl  _  ^1353  -36  57  57 

57  _  a/I353  +  21  _  q  a/I353  —  27 

X2  _  ^1353  -  21  _  16  +  16 


(3.3) 


The  second  fraction  in  each  step  is  found  by  rationalizing.  At  each  step,  the  integers 
taken  out  are  bi  and  the  remaining  fractions  are  between  0  and  1.  After  subtracting  bt  the 
remaining  fraction  is  inverted  to  find  xl+\ .  As  a  point  of  reference,  we  have  approximated 
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so  far  that  a/1353  ~  36  +  — ^j.  SQUFOF  stops  here  because  16  is  a  perfect  square.  Taking 
the  conjugate  of  the  numerator  and  the  square  root  of  the  denominator,  we  obtain: 


,  ^1353  +  27  a/I353  -  33 

x>= - 4 - =  15+ - 4 - 

,  4  a/I353  +  33  ,  \/l353  —  33 

x  =  ,  - = - =  H - (3.4) 

a/1353  -  33  66  66 

Here,  since  the  residue  33  is  repeated,  we  quickly  find  that  33  is  a  factor  of  1353.  1353  = 
33-41  =  3-11-41. 

Shanks  developed  Square  Forms  Factorization,  or  SQUFOF,  based  on  a  concept  called 
the  infrastructure  of  the  class  group.  This  refers  to  the  relationship  between  distance  and 
composition.  Define  the  distance  from  Frn  to  Fn  by: 


D(Fm,  Fn)  =  log(  xk)  (3.5) 

k=m-\- 1 

where  xk  are  the  corresponding  terms  in  the  continued  fraction  expansion.  Since  distance 
is  the  log  of  the  product  of  the  terms  in  between,  roughly  speaking,  distance  is  nearly 
proportional  to  the  difference  in  indices.  Then  the  location  of  the  form  resulting  from  the 
composition  of  two  quadratic  forms  is  controlled  by: 


D(F1  *  GuFn  *  Gm)  =  D(Fi,  Fn)  +  D{GU  Gn)  +  C  (3.6) 

where  Q  is  small3.  This  formula  is  proven  in  the  arguments  leading  up  to  Theorem  A. 5. 2. 

By  equation  (3.6),  the  distance  from  a  symmetry  point  is  doubled  when  a  quadratic  form 

is  squared,  so  that  the  square  root  of  a  quadratic  form  is  half  the  distance  from  a  symmetry 

point  as  the  original  quadratic  form.  The  change  made  to  the  continued  fraction  expansion 

3The  bound  for  £  is  proportional  to  log  log  N  while  the  two  distances  on  the  right  side  of  the  formula 
have  a  bound  proportional  to  log  N ,  so  (  is  usually  negligible,  although  it  may  be  calculated  if  necessary. 
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upon  the  occurrence  of  a  perfect  square  corresponds  exactly  to  taking  the  square  root  of 
a  quadratic  form  and  reversing  its  direction,  so  that  in  the  second  phase  (after  the  change 
in  sequence  is  made)  the  algorithm  is  going  backwards  on  a  different  cycle.  The  distance 
covered  from  making  this  change  until  obtaining  the  symmetry  point  will  be  one  half  the 
distance  covered  before  finding  this  square.  Since  this  distance  may  be  known  as  well  as 
necessary,  it  is  possibly  to  use  some  of  the  quadratic  forms  found  along  the  way  to  shorten 
this  return. 

Formally,  here  is  the  version  of  Square  Forms  Factorization  that  Shanks  intended: 

Given  N  to  be  factored: 

Qo<-l,P0<-  \_VN\1Ql^N-P$ 

r  <-  lVN\ 

while  Qt  7^  perfect  square  for  some  %  even 


L  J 

Pi  <-  hQi  ~  Pi-i 

Qi+l  Qi-1  +  bi(Pi- 1  —  Pi) 

if  i  —  2n  for  some  n 

_ _  D2  _  AJ 

Store  (Qi,  2  •  P/)  F0  =  (VQu  2  ■  P-i,  ) 

Compose  Fq  with  stored  forms  according  to  the 
binary  representation  of  i/2  and  store  result  to  Fn. 

F0  =  (A,B,C ) 

Q0  \A\,P0  <-  B/2,Q1  \C\ 

Qo  <-  Qi,Po  <-  Po,qi  <-  Qo 

while  Pi  ^  Pi- 1  and  pi  ^  p%-\ 

Apply  same  recursive  formulas  to  (<3tu-Po>Qi)  an(l  (qo,Po,Qi) 
If  Pi  =  Pi_ i,  either  Qi  or  Qi/ 2  is  a  nontrivial  factor  of  N. 

If  pi  =  Pi- 1,  either  qi  or  qi/ 2  is  a  nontrivial  factor  of  At. 


Figure  3.2:  SQUFOF  with  Composition 


Note  that  it  is  faster  to  approximate  roughly  which  quadratic  forms  are  needed  to  find 
the  symmetry  point  than  to  actually  calculate  the  distance  exactly.  The  result  is  that  after 
this  composition,  the  quadratic  form  is  close  to  the  symmetry  point  but  not  exactly  on  it. 
In  the  last  while  loop,  the  search  for  the  symmetry  point  has  to  be  made  in  both  directions, 
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but  this  search  is  typically  very  short.  The  only  forms  that  need  to  be  stored  are  those  with 
an  index  that  is  a  power  of  two.  The  predictions  of  which  squares  provide  a  non-trivial  factor 
is  unnecessary  since  with  Fast  Return  it  is  faster  to  test  each  square  found. 

For  example,  if  a  square  were  found  on  step  56,  there  should  be  roughly  28  steps  backward 
in  the  second  phase.  If  after  taking  the  square  root  and  reversing  the  direction,  Fq  is 
composed  with  the  quadratic  form  with  index  16,  the  form  with  index  8,  and  the  form  with 
index  4,  the  result  will  be  very  close  to  the  symmetry  point. 

3.3  After  Shanks 

Shanks  died  in  1996  leaving  several  papers  on  Square  Forms  Factorization  in  unpublished 
form  [27].  These  were  incomplete  and  contained  very  few  proofs.  However,  Shanks  suggests 
he  had  proofs.  In  his  main  paper,  for  example,  he  states  that  SQUFOF  has  an  average 
runtime  of  0(yfN)  and  says  this  will  be  proven  later.  Since  the  paper  was  not  complete,  it 
may  never  be  known  whether  or  not  Shanks  possessed  a  proof  of  this  critical  detail. 

Although  Shanks  did  not  publish  his  work,  it  has  been  picked  up  by  several  people. 
Williams  [26],  in  1985,  provides  a  description  of  SQUFOF  based  on  ideals  and  lattices. 
Williams  demonstrates  how  continued  fractions  can  be  seen  as  an  intuitive  extension  of  the 
theory  of  lattice  minima  (§4)  and  uses  lattices  to  connect  composition  with  distance  (§6). 
Williams  went  further  by  generalizing  these  results  complex  cubic  fields.  Riesel  [18],  also  in 
1985,  gives  a  good  exposition  of  SQUFOF  using  the  theory  of  continued  fractions.  Buell  [2], 
in  1989,  provides  a  thorough  analysis  of  quadratic  forms  and  gives  a  description  of  SQUFOF. 

Recently,  the  thesis  of  Jason  Gower  in  [9]  discussed  the  runtime  analysis  of  SQUFOF  and 
Shanks’  belief  that  the  expected  runtime  is  0(\fN)  =  L(l,  1/4).  Jason  Gower  provides  a 
good  description  of  continued  fractions,  quadratic  forms,  and  the  basic  SQUFOF  algorithm. 
He  also  did  some  research  and  experiments  on  the  use  of  a  multiplier  and  discovered  that 
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it  could  be  an  effective  method  to  accelerate  the  algorithm.  His  analysis  of  the  runtimes  is 
based  on  the  distribution  of  square  forms  through  the  ambiguous  cycles  and  is  supported  by 
extensive  data.  This  work  was  completed  in  December  of  2004. 

However,  although  much  of  the  theory  about  SQUFOF  was  known,  some  of  it  was  not 
and  most  of  it  was  very  scattered  with  few  proofs.  Until  now,  no  complete  explanation  of  why 
SQUFOF  works  had  been  published.  Therefore,  in  the  appendices  I  have  compiled  all  of  this 
information  and  provided  proofs  and  examples.  I  first  analyzed  continued  fractions,  showing 
that  taking  the  conjugate  of  the  numerator  provide  a  new  and  intuitive  explanation  of  the 
reversal  of  direction  (Lemmas  A. 1.3  -  A. 1.4).  This  understanding  then  provides  immediate 
results  about  the  periodicity  and  symmetries  of  the  cycles  (Lemmas  A. 1.5  -  A. 1.7).  These 
symmetry  points  correspond  to  symmetry  points  in  the  quadratic  form  cycle  (Lemma  A. 2. 6). 

The  fact  that  the  square  of  any  ambiguous  cycle  is  the  principal  cycle  is  of  fundamental 
importance  to  SQUFOF  and  is  proved  in  Theorem  A. 2. 8.  This  result  was  probably  under¬ 
stood  by  Shanks,  but  apparently  never  stated  explicitly.  This  fact  implies  that  after  taking 
the  square  root  of  a  quadratic  form  in  the  principal  cycle  the  result  is  in  an  ambiguous  cycle 
so  that  there  will  be  a  symmetry  point. 

The  application  of  ideals  to  quadratic  forms  required  a  slight  generalization  to  the  stan¬ 
dard  description  of  ideals  (Lemma  A. 3. 2).  This  slight  extra  generality  allows  a  proof  that 
multiplication  of  ideals  corresponds  to  composition  of  quadratic  forms.  The  well-known  re¬ 
sult  about  the  multiplication  of  ideals  was  also  generalized  slightly  (A. 3. 3)  so  that  the  case 
of  quadratic  forms  with  discriminant  =  1  (mod  4)  could  be  handled  without  working  in 
Z  y/^+1  .  The  connection  between  lattice  reduction  and  continued  fraction  reduction  also 
shows  that  a  distance  for  reduction  may  be  computed  (Lemma  A. 5.1). 

Theorem  A. 5. 2  relates  composition  of  quadratic  forms  with  distance  and  is  a  generaliza¬ 
tion  of  the  distance  formula  provided  by  Williams  in  [26].  This  generalization  is  essential  for 
SQUFOF  as  it  relates  distances  in  different  non-principal  cycles  with  composition.  Lemma 
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A. 5. 3,  which  relates  the  entire  distance  around  different  cycles,  is  proved  based  on  this.  With 
these  facts,  it  is  possible  to  prove  that  the  distance  traveled  in  the  non-principal  cycle  before 
finding  a  symmetry  point  will  be  one  half  of  the  distance  traveled  in  the  principal  cycle  to 
find  the  square  form. 

Appendix  A. 6  provides  a  complete  description  of  Square  Forms  Factorization,  including 
an  explanation  of  why  Shanks  was  able  to  know  whether  or  not  a  square  was  proper  and  an 
explanation  of  what  Shanks  referred  to  as  Fast  Return.  Fast  Return,  an  idea  that  is  typically 
not  implemented  with  SQUFOF,  uses  composition  to  significantly  speed  up  the  second  phase 
of  the  algorithm. 
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Chapter  4 


Computer  Based  Experiments 


Appendix  B  includes  the  C  implementation  of  SQUFOF  with  Fast  Return.  These  algo¬ 
rithms  have  also  been  implemented  in  Magma.  These  have  both  been  tested  extensively 
for  functionality.  The  runtimes  have  also  been  compared  to  elliptic  curve  method,  multiple 
polynomial  quadratic  sieve,  and  the  number  held  sieve.  Although  a  thorough  comparison  of 
SQUFOF  with  other  algorithms  has  not  been  done,  it  is  fairly  clear  that  on  a  single  computer 
SQUFOF  is  slower  than  elliptic  curves,  the  number  held  sieve,  or  the  multiple  polynomial 
quadratic  sieve  for  most  numbers.  However,  there  have  been  several  other  areas  that  have 
been  researched.  It  has  been  conjectured  by  Pomerance  [17]  that  a  parallel  implementation 
of  SQUFOF  would  be  competitive.  This  area  of  research  provided  some  valuable  insights 
into  these  algorithms. 

4.1  Two  Parallel  Implementations 

The  current  trend  in  ultra-fast  computing  is  to  use  a  large  number  of  inexpensive  processors 
instead  of  a  single  expensive  processor  [7].  With  the  large  amount  of  computation  required 
for  factorization,  the  efficiency  of  a  parallel  implementation  is  especially  important  for  these 
algorithms.  In  1982,  Shanks  and  Cohen  attempted  a  parallelization  by  having  multiple  pro- 
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cessors  attempt  to  factor  N,  3N,  5 N,  etc.  Gower’s  recent  Ph.D.  thesis  (under  Dr.  Wagstaff) 
[9]  analyzed  the  use  of  multipliers  and  found  them  to  be  effective  in  general  but  didn’t  provide 
much  evidence  on  their  efficiency  for  parallelization. 

Parallelization  by  multipliers  has  been  implemented  in  C  using  the  MPI  library  for  the 
purpose  of  runtime  comparison  and  is  included  in  Appendix  B. 

This  research  has  discovered  that  composition  of  quadratic  forms  provide  a  different 
possible  approach  to  parallelization.  By  computing  a  quadratic  form  several  steps  into  the 
cycle  and  squaring  it  a  reasonable1  number  of  times,  a  quadratic  form  far  out  into  the  cycle 
can  be  found.  Call  it  F.  The  first  processor  may  be  assigned  to  search  this  segment  for  a 
perfect  square.  The  next  processor  can  search  from  F  to  F2,  the  next  from  F2  to  F3,  etc. 
Each  of  these  segments  will  be  roughly  equal  sizes.  When  a  processor  finds  a  perfect  square, 
it  may  use  the  quadratic  forms  involved  in  computing  its  segment  in  order  to  come  close  to 
the  symmetry  point.  Specifically,  this  last  step  is  made  easier  if  the  starting  point  of  the 
segment  is  found  by  squaring.  In  pseudo-code: 

Given  N  to  be  factored: 
r  <-  [VN\ 

F0  <—  (1, 2 r,  N  —  r 2) 

Cycle  F0  several  steps  forward. 

for  i  —  1  to  size  (size  is  the  logarithmic  size  of  a  segment.) 

Fi  Fi_i  *  Fi_ i 

Processor  0: 

Assign  one  processor  to  search  from  Fq  to  Fslze . 


F start 

-  F 

±  size 

Fend  * 

F2 

size 

FrootS  < 

—  F 

±  size— 

FrootE 

( —  F  • 

*  L  size 

F step 

F size—  1 

while  A  factor  hasn’t  been  found 

Wait  for  a  processor  to  be  free  and  send  Fstart ,  Fendl  and  Froots. 
Fstart  *  Fend 
FrootS  ‘  Frooi  jr 


Empirically,  20-30  times  works  well.  This  parameter  isn’t  critical. 
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F rootE  *  F rootE  *  F sfep 

rp  _  t?2 

1  end  1  rootE 

Processor  n : 

Receive  FS£art ,  Fenc[ ,  &nd  Frootg 
count  0 

F0  =  (A,B,C) 

while  A  factor  is  not  found  and  Fstart  ^  Fend, 

Cycle  F start  forward  2  steps, 
count  <—  count+1 
if  A  is  a  perfect  square 

-p  , _  p-1/2 

rtest  r  start 

Ftest  *  FfCst  *  FrootS 
for  j  =  size  to  1 
if  count  >  2j 
Ftest  *  Ftest  *  Fj 
count  count  —  2-7 
Search  in  both  directions  from  Ftest  for  a  symmetry  point, 
if  Factorization  found  at  symmetry  point,  output  and  quit, 
if  A  factor  is  still  not  found 

Receive  new  Fstart,  Fend ,  and  Froots  and  start  over. 


(This  loop  composes  Ftest  with  the  necessary 
forms  to  bring  it  close  to  the  symmetry  point.) 


Figure  4.1:  Parallel  SQUFOF 


Since  there  is  no  overlap  between  the  segments  searched  by  the  processors  and  since 
the  perfect  squares  are,  apparently,  distributed  randomly,  this  parallelization  should  be  ex¬ 
tremely  efficient.  The  size  of  the  segments  is  determined  by  how  far  the  first  for  loop 
goes.  There  are  only  two  hazards  when  choosing  this.  If  the  segment  size  is  too  small,  the 
processors  will  finish  their  segments  so  quickly  that  receiving  new  segments  will  become  a 
bottleneck.  Alternately,  if  the  segments  are  too  long,  the  processors  may  divide  up  more 
than  the  entire  cycle,  so  that  there  is  overlap.  However,  except  for  rare  numbers  that  will 
factor  trivially  fast  regardless,  there  is  significant  room  in  between  these  two  bounds. 

This  parallel  version  has  been  implemented  in  C  using  the  MPI  library.  The  code  is 


included  in  Appendix  B. 
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There  are  several  characteristics  of  these  algorithms  that  were  analyzed  in  this  compar¬ 
ison.  The  test  integers  were  all  products  of  randomly  chosen  primes  of  roughly  equal  size 
generated  by  Magma.  80  bit,  100  bit  and  120  bit  integers  were  all  tested  on  20,  30,  40  and 
50  processors.  This  allows  an  analysis  of  both  how  each  algorithm  is  affected  by  the  size  of 
the  integers  and  how  efficiently  each  algorithm  uses  an  increasing  number  of  processors. 

The  first  comparison  is  an  analysis  of  which  algorithm  was  faster  on  a  larger  number  of 
the  integers.  Table  4.1  shows  these  results. 


Processors 


20 

30 

40 

50 

80 

0.39 

0.38 

0.39 

0.30 

bits 

100 

0.37 

0.34 

0.43 

0.40 

120 

0.35 

0.31 

0.29 

0.39 

Table  4.1:  %  Wins  by  the  Segments  Parallelization 


The  data  is  clearly  showing  that  parallelization  with  multipliers  is  faster  than  segments 
and  shows  no  apparent  pattern  as  to  how  the  size  of  the  integers  or  the  number  of  the 
processors  might  affect  this. 

The  next  comparison  is  an  analysis  of  the  average  runtimes.  This  comparison,  demon¬ 
strated  in  Table  4.2  is  a  bit  more  interesting. 


Bits 


Segments _  _ Multipliers 


Processors: 

80 

100 

120 

80 

100 

120 

20 

5.04 

157.27 

5416.97 

2.95 

82.62 

3255.65 

30 

3.03 

116.76 

4287.37 

1.95 

59.75 

2028.72 

40 

2.53 

74.50 

2690.76 

1.68 

52.90 

1683.57 

50 

2.13 

77.67 

2538.80 

1.23 

53.60 

1287.80 

Table  4.2:  Comparison  of  Average  Runtimes 


A  quick  survey  of  this  data  shows  that  for  the  segments  parallelization,  the  average 
runtime  was  cut  in  half  from  20  processors  to  40,  while  the  multipliers  implementation 


38 


didn’t  do  quite  so  well.  This  can  be  represented  by  multiplying  each  runtime  by  the  number 


of  processors  and  then  dividing  by  the  mean  runtime  for  that  many  bits  to  place  all  the 


results  on  a  single  graph: 


Segments 


-•-80  Bit 
^100  Bit 
-*-120  Bit 
-x- Average 


#  Processors 


Figure  4.2:  Segments  Efficiency  of  Parallelization 


Figure  4.3:  Multipliers  Efficiency  of  Parallelization 


These  graphs  show  that  the  efficient  use  of  multiple  processors  for  the  segments  paral¬ 
lelization  is  roughly  unaffected  by  increasing  the  number  of  processors,  while  the  multipliers 
parallelization  is  less  efficient  at  using  a  larger  number  of  processors.  This  is  the  expected 
result.  As  Jason  Gower  demonstrated  in  [9],  the  use  of  a  multiplier  can  decrease  the  average 
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runtime  by  an  average  of  27%.  Therefore,  for  small  numbers  of  processors,  using  multi¬ 
pliers  should  immediately  cut  the  average  runtime  down.  However,  for  larger  numbers  of 
processors,  the  multipliers  available  aren’t  used  as  efficiently. 

This  same  trend  can  be  expressed  by  graphing  the  ratios  of  average  runtimes: 


Figure  4.4:  Ratios  of  Segments  v.  Multipliers 

Although  the  data  isn’t  completely  clear,  the  trend  is  toward  segments  being  faster  than 
multipliers  if  enough  processors  are  used.  Based  on  the  ratios,  a  linear  regression  predicts 
a  crossover  at  180  processors.  However,  as  larger  multipliers  are  used,  the  efficiency  of  the 
multipliers  should  drop  off  significantly,  so  the  actual  crossover  will  be  earlier  than  that. 

Considering  that  directly,  SQUFOF  probably  has  average  runtime  of  0{\fN ),  and  it 
appears  to  use  the  processors  available  with  near  perfect  efficiency,  the  expected  runtime 
on  multiple  processors  is  0(\/rN/np),  where  np  is  the  number  of  processors.  Based  on  the 
data  that  for  120  bits,  the  factorization  took  an  average  of  31.91  processor  hours,  it  is 
then  expected  that  RSA-576,  the  last  RSA  challenge  number  factored,  would  take  31.91  * 

2(576— 120)/4  =  g  gg  *  1035 

processor  hours.  Thus,  for  example,  with  1020  processors,  it  would 
take  6.63  *  1015  hours  =  7.57  *  1011  years.  Both  of  these  figures  are  extremely  unfeasible. 
Therefore,  a  parallel  implementation  of  SQUFOF  will  probably  never  be  competitive  with 
other  algorithms. 
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4.2  Test  of  Direction 

Based  on  the  symmetry  of  the  cycles  of  quadratic  forms,  any  fast  test  to  determine  whether 
or  not  the  continued  fraction  expansion  is  in  the  correct  direction,  that  is,  whether  or  not  it 
has  passed  the  factorization  yet,  would  provide  a  faster  factorization  algorithm  by  performing 
a  binary  search  (Figure  4.5). 


Given  N  to  be  factored: 

Qo  <-  1,  Po  <-  IV^VJ 

Apply  equation  (2.2)  for  2  steps. 

Fo  ( Q2 ,  2  •  P2,  —Q3) 

i  <-  0 

while  Fi  is  in  the  right  direction 

Fi+ 1  <-  Fi*  Fi 

i  *  +  1 

Flast  *  F{—] 

while  i  >  0 

if  Fiast.  *  Fi  is  in  the  right  direction 
Flast,  *  Fiast  *  Fi 
i  <—  i  —  1 

Search  forward  from  Fiast  (should  be  only  a  couple  steps) 
to  obtain  the  factorization. 


Figure  4.5:  Binary  Search  based  on  a  Test  of  Direction 


Exatnple  8.  N  =  2035153,  Q0  —  1,  P0  —  1426 

Then,  Fo  =  (1,  2  •  1426,  —1677).  The  square  of  this  form  would  be  itself,  so  we  recursively 
step  forward  a  few  steps  to  obtain:  F$  =  (663,  2774,  —168). 

This  form  is  squared  until  it  is  past  the  symmetry  point: 

F8*F8  =  (1569, 1522,  -928)  =  F22 
F22  *  F22  =  (896,  2798,  -87)  =  F44 
Fm  *  Fm  =  (1648, 1726,  -783)  =  F82 
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F82  *  F8 2  =  (411,  2104,  —2259)  is  past  the  symmetry  point. 

F82  *  Fu  =  (9,  2846,  -1136)  =  F134 
F\ 34  *  F22  is  past  the  symmetry  point. 

F134  *F8  =  (1153, 1898,  -984)  =  FU4 

Recursively  stepping  forward  from  here,  it  is  quickly  found  that  F147  =  (—1008,  20  1  8, 1009) 
and  Fi48  =  (1009,2018,-1008).  Therefore,  1009  is  the  symmetry  point  and  2035153  = 
1009  •  2017. 

The  decisions  for  this  example  of  whether  or  not  a  form  was  reversed  (past  the  symmetry 
point)  were  determined  by  merely  comparing  with  the  actual  continued  fraction  expansion. 
However,  ideally  this  could  be  done  without  doing  the  entire  expansion.  The  function  that 
determines  whether  or  not  a  form  is  reversed  is  called  a  test  of  direction.  There  were  several 
possible  candidates: 

Conjecture  1.  If  QfQi-i,  (Qi)3  \  Qi-i,  and  Qi  is  not  a  power  of  2,  then  the  continued  fraction 
expansion  is  in  the  correct  direction. 

In  all  other  cases,  this  test  provides  no  information.  The  condition  Qi\Qi-\  does  not 
occur  very  frequently,  so  on  its  own  it  cannot  provide  a  useful  test  of  direction.  However, 
the  attempt  to  explain  this  empirical  pattern  has  provided  several  possible  alternatives: 

First,  it  is  possible  to  find  a  multiple  k  such  that  in  the  expansion  with  Qi  unchanged, 
Pi  replaced  by  kPi  and  N  replaced  by  k'2N,  this  condition  is  met.  For  notation,  let  Qi  = 
Q Q.  Since  using  —k  instead  of  k  produces  the  same  sequence  in  the  opposite  direction, 
it  is  necessary  to  relate  this  sequence  to  the  original  sequence  in  order  to  extract  useful 
information.  From  this,  it  has  so  far  been  found  that  if  Qi+3|Qi  and  Q'0\Q'_i,  then  the  original 
expansion  is  in  the  correct  direction.  Finding  multiples  such  that  the  second  condition  is 
satisfied  can  be  accomplished  by  using  continued  fractions  to  find  the  convergents  of  , 

where  Ri  =  yfN  =  Pi  —  PjpN  (mod  Qj). 
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This  provides  a  test  of  direction  that  returns  a  result  about  2%  of  the  time.  This  would 
potentially  provide  a  polynomial  time  factorization  algorithm,  except  that  the  time  required 
to  do  this  test  of  direction  increases  with  the  size  of  N. 

One  other  possible  test  of  direction  involves  reduction  distance.  After  a  quadratic  form 
is  squared,  it  often  requires  multiple  steps  in  order  to  reduce  to  a  quadratic  form  within 
the  standard  bounds.  This  is  referred  to  as  reduction  distance.  Since  the  overall  distance 
from  a  symmetry  point  is  doubled  when  a  quadratic  form  is  squared,  it  makes  some  intuitive 
sense  that  this  distance  from  the  nearest  symmetry  point  should  be  related  to  this  reduction 
distance.  Empirically,  the  reduction  distance  for  a  quadratic  form  oriented  away  from  the 
nearest  symmetry  point  tends  to  be  shorter.  Concerning  the  conjecture,  if  Qt  \  Qi-i,  then 
the  square  of  (Qi,  Pi-i,  —  Qi-i)  is  (Q2,  Pi-i,  —  Qi-i/Qi)  and  the  reduction  distance  is  0,  so 
this  would  provide  an  explanation  for  Conjecture  1. 

Empirically,  comparing  reduction  distances  is  a  test  of  direction  that  generally  works  near 
the  symmetry  point,  but  is  inconclusive  away  from  the  symmetry  point.  Thus,  one  other  idea 
may  be  necessary  in  order  to  apply  this.  Using  the  class  number  formula  [3],  it  is  possible  to 
approximate  what  the  distance  to  the  symmetry  point  should  be.  This  approximation  may 
be  able  to  come  close  enough  to  the  symmetry  point  that  some  combination  of  the  two  tests 
of  direction  may  be  able  to  provide  the  necessary  information. 

4.3  Fast  Numbers 

The  other  area  of  research  was  whether  there  were  some  integers  that  would  factor  extremely 
fast  and  whether  these  would  pose  a  security  risk.  In  general,  it  seems  likely  that  the  numbers 
that  factor  extremely  fast  are  of  the  form  ( f2m  +  c)(g2m  +  d),  where  m  is  “large”  and  a,b,c, 
and  d  are  relatively  small.  The  general  conditions  on  these  variables  appear  to  be  very 
complicated.  The  simpler  special  case  of  p(p  +  2a),  for  which  it  is  fairly  easy  to  analyze  the 
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conditions  exactly  and  determine  how  often  this  case  occurs,  provides  some  indication  of  the 
complexity  of  this  question.  Although  an  exact  solution  may  be  developed  for  any  particular 
case,  for  the  general  case  it  has  so  far  only  been  possible  to  approximate  how  often  these 
cases  occur. 

Proposition  1.  Let  N  =  p(p  +  2a),  where  p  and  a  are  integers  with  0  <  a  <  v/p.  Then  the 
Square  Forms  Factorization  algorithm  will  provide  a  factorization  of  N  on  the  second  step 
with  runtime  O(logiV).  Specifically,  Q2  =  a2  and  Q\  ^  a. 

Proof:  Let  Qo  —  1. 

(a  —  l)2  <  a2  <  p  <  2p 
(a  -  l)2  -  2p  <  0 

p2  +  2 ap  +  (a  —  l)2  —  2p  <  p2  +  2 ap  <  p2  +  2 ap  +  a2 
{p  +  a  —  l)2  <  N  <  {p  +  a)2 
iViVj  =  p+  a  -  1 

Therefore,  Pq  =  p  +  a  —  1  and  Q\  =  N  —  P02  =  2p  —  (a  —  l)2.  b\  =  •  Now, 

a2  —  a  —  p  <  0 
p  +  a  —  1  <  2p  —  (a  —  l)2 

2p+2a— 2  9 

2p— (a— l)2  ^  Z 

bi  =  1 

Therefore,  P1  =  p  ~  a(a  -  1).  Q2  =  =  2a~2ppZa{a-ipy  =  °2- 

If  Q i  =  a,  then  2p  —  (a  —  l)2  =  a,  so  that  a  =  1±v/|p~3  and  8p  —  3  =  5  (mod  8)  and  is 
thus  not  a  perfect  square.  Since  a  is  an  integer,  this  is  then  impossible. 

VN  -  P  +  a(a  -  1)  a fN  —  p 

Vo  —  —  («  ~  1)  H 

a  a 
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a  VN  +  p  y/N  —  p 

Vl  A /N  -  p  2 p  +  2 p 

The  factor  of  p  is  thus  found  extremely  rapidly.  Since  the  only  part  of  the  algorithm  that 
takes  a  significant  amount  of  time  is  approximating  the  square  root  of  N,  the  algorithm  has 
a  runtime  of  0(log  N)  for  numbers  that  meet  this  criteria.  QED 

For  the  other  cases,  the  fact  that  p(4p  +  4a  +  1)  and  p(Ap  +  4a  +  3)  have  significantly 
different  bounds  for  a  should  be  some  indication  as  to  how  quickly  this  problem  can  get 
messy,  ffowever,  some  approximations  can  still  be  determined  empirically.  Let  N  =  ( f2m  + 
c)(g2m  +  d).  Then  empirically  N  factors  quickly  for  approximately  the  same  cases  as  when 
f2g2N  =  ( f2g2m  +  g2c )  *  ( f2g2m  +  f2d )  factors  quickly.  Let  p  =  f2g2m  +  g2c )  and  let 
a  =  (f2d  —  g2c)/ 2  or  (f2d  —  g2c—  l)/2,  then  f2g2N  =  p(p  +  2a)  or  p(p  +  2a  +  1,  respectively. 
In  either  case,  the  two  factors  (f2g2m  +  g2c)  and  (f2g2m  +  f2d)  must  have  roughly  half  their 
bits  the  same. 

For  a  given  N  =  qr,  consider  the  probability  that  there  exist  /  and  g  such  that  g2q  and 
f2r  have  the  same  number  of  bits  and  the  first  3  bits  are  the  same.  This  probability  can 
easily  be  tested  empirically.  The  continued  fraction  expansion  for  q/r  provides  potential 
values  for  /  and  g.  A  test  of  1015  pairs  of  100  bit  integers  and  1015  pairs  of  300  bit  primes 
show  that  this  happens  approximately  11%  of  the  time  in  each  case,  with  only  a  couple  of 
pairs  that  had  multiple  values  of  /  and  g  meeting  the  criteria.  It  seems  reasonable  to  assume 
that  if  the  size  of  the  paired  integers  wasn’t  required  to  be  equal,  the  probability  would  be 
different  but  would  still  be  fairly  constant.  Let  this  constant  be  k.  Then  for  each  bit  after 
that,  there  is  a  50%  probability  of  q  and  r  having  the  same  value.  Since  roughly  the  first 
half  of  the  bits  must  be  equal,  the  probability  is  k  *  (,5)log2 N/2~3  =  8 k/y/N  =  0(1/VN). 
Thus,  these  numbers,  although  interesting,  are  extremely  rare  and  should  not  be  a  concern 
in  normal  operations. 
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Chapter  5 


Conclusions  and  Future  Work 


This  report  has  described  and  explained  why  Square  Forms  Factorization  works,  completing 
the  work  started  by  Daniel  Shanks.  Specifically,  a  thorough  correspondence  was  described 
between  quadratic  forms,  continued  fractions,  ideals  and  lattices,  resulting  in  a  proof  of  a 
generalized  distance  formula.  This  allowed  an  explanation  of  how  Square  Forms  Factoriza¬ 
tion  works.  Several  variations  have  been  considered.  An  analysis  shows  that  the  use  of 
multipliers  decreases  the  average  runtime  by  25%.  However,  the  use  of  multipliers  alone 
is  an  inefficient  method  of  implementing  Square  Forms  Factorization  on  a  large  number  of 
processors.  Parallelization  by  dividing  the  cycle  into  segments  for  each  processor  appears 
to  be  more  efficient  with  larger  number  of  processors.  The  crossover  between  these  two 
methods  appears  to  be  at  about  100  processors.  In  light  of  this,  the  best  parallelization 
would  probably  involve  some  combination.  Two  alternatives  were  considered  but  not  tried. 
It  would  be  fairly  simple  to  use  a  single  multiplier  and  then  use  segments  to  parallelize.  It 
may  be  even  more  efficient  to  use  several  multipliers  and  then  use  segments  to  divide  each 
of  those  between  many  processors. 

With  the  parallelization,  it  is  important  to  realize  that  these  parallelizations  are  all  very 
easy  to  implement  and  would  be  very  effective  as  a  distributive  process.  The  only  information 
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that  the  master  processor  needs  to  send  is  basic  information  needed  to  assign  a  segment. 
The  other  processor  then  works  independently  until  it  either  finds  a  factor  or  finishes  its 
segment. 

Another  alternative  considered  using  a  test  of  direction  to  perform  a  binary  search.  How¬ 
ever,  the  test  of  direction  only  appeared  to  be  effective  near  the  symmetry  points  and  even 
with  an  approximation  of  the  regulator  it  was  not  possible  to  get  close  enough  to  the  sym¬ 
metry  points  for  the  test  of  direction  to  be  useful. 

One  other  property  of  Square  Forms  Factorization  that  was  considered  is  which  numbers 
factor  extremely  fast.  Although  there  are  integers  of  every  size  that  factor  in  time  0(log  N), 
they  are  rare.  The  result  was  that  an  integer  N  has  a  probability  0{l/y/N)  of  meeting  this 
criteria.  Therefore,  this  isn’t  a  threat  worth  considering  seriously. 
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Chapter  6 


Glossary 


Conjugate:  modification  to  an  algebraic  number  changes  the  sign  of  a  certain  term  in  an 
expression,  traditionally  the  part  of  an  expression  that  is  either  irrational  or  imaginary. 
In  the  context  of  continued  fractions,  it  applies  to  changing  the  sign  of  a  residue,  the 
integer  part  of  the  numerator  of  an  x,  in  the  continued  fraction  expansion. 

Cryptography:  the  process  of  encoding  and  decoding  information  to  prevent  it  from  being 
read  by  anyone  who  intercepts  the  message.  Used  in  variety  of  civil  and  military 
applications. 

Euclid’s  Algorithm:  a  fast  algorithm  for  determining  the  greatest  common  divisor  of  two 
integers.  Given  x  and  y ,  the  extended  Euclidian  algorithm  also  determines  the  coeffi¬ 
cients  a  and  b  such  that  ax  +  by  =  gcd(x,  y). 

Floor:  the  greatest  integer  less  than  or  equal  to  a  given  number.  Symbol:  [x\ 

Greatest  Common  Divisor:  largest  integer  that  divides  a  pair  or  group  of  integers.  Sym¬ 
bol:  gcd(x,y ) 

Modular  Arithmetic:  two  integers  are  considered  equal  (congruent)  if  their  difference  is 
divisible  by  a  base.  Thus,  3  =  10  (mod  7).  Numbers  are  represented  by  integers 
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between  0  and  N  —  1,  where  N  is  the  base.  Multiplication,  addition,  and  subtraction 
are  normal,  except  that  the  results  are  reduced  to  the  range  0  to  N  —  1.  Division  is 
performed  by  reducing  fractions  to  least  terms,  applying  an  extended  Euclid’s  algorithm 
to  find  the  inverse  of  the  denominator,  and  then  performing  multiplication.  Symbol: 
(mod  N) 

Modulo:  operation  related  to  division  that  returns  the  remainder: 
yj  =  6jy,  so  73  modulo  11  =  7.  Symbol:  %  or  (mod  N). 

Non-square  integer  that  is  not  a  perfect  square  of  another  integer. 

NUCOMP-NUDUPL:  algorithms  designed  by  Daniel  Shanks  to  perform  composition  of 
quadratic  forms  quickly. 

Parallelization:  the  distribution  of  a  single  task  between  multiple  processors  working  si¬ 
multaneously. 

Perfect  Square:  integer  that  is  the  square  of  another  integer.  Thus,  9  is  a  perfect  square 
because  32  =  9. 

Prime:  integer  that  has  no  divisors  other  than  itself  and  1.  The  first  several  are  2,  3,  5,  7, 11.... 

Primitive:  a  characteristic  that  implies  not  having  common  factors.  For  quadratic  forms, 
a  form  whose  coefficients  do  not  contain  a  common  factor.  For  ideals,  an  ideal  whose 
elements  does  not  contain  a  common  integer  factor. 

Pseudo-square:  the  denominator  of  an  xt  in  the  continued  fraction  expansion,  denoted  (),. 
When  Q0  =  1,  (— 1  )lQj  is  a  quadratic  residue  and  in  general  —  QiQi-i  is  a  quadratic 
residue. 

Quadratic  Reciprocity  Law  (Gauss):  determines  which  numbers  are  quadratic  residues 


of  a  prime: 
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Symbol:  =  1  if  x2  =  a  (mod  p)  has  a  solution,  —1  if  it  does  not,  and  0  if  p  \  a. 

Theorem:  For  p  and  q  distinct  primes: 

\qj  p 

(p)  =  (-i)«',2-i)/8 

(y)  =  (-D<”-I,/2 

Quadratic  Residue:  a  perfect  square  modulo  some  base  N.  2  is  a  quadratic  residue  of  7 
because  32  =  2  (mod  7). 

Relatively  Prime:  a  description  of  two  integers  that  have  no  common  divisors.  Thus,  8 
and  15  are  relatively  prime,  even  though  they  are  not  prime  by  the  normal  definition. 

Residue:  an  integer  that  remains  in  the  numerator  after  bt  has  been  subtracted  from  ay  in 
the  continued  fraction  expansion. 

RSA:  a  cryptology  algorithm  named  after  Rivest,  Shamir,  and  Adleman.  It  was  earlier 
developed  by  Clifford  Cooks  of  GCHQ,  but  this  was  only  recently  declassified.  Its 
security  is  dependent  on  the  difficulty  of  factorization  [5]. 

SQUFOF:  Square  Forms  Factorization  algorithm,  developed  by  Daniel  Shanks  in  1975. 
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Appendix  A 
Proofs 


The  largest  part  of  this  research  has  been  developing  formal  proofs  of  the  concepts  sup¬ 
porting  Square  Forms  Factorization.  Section  A.l  introduces  and  proves  a  number  of  results 
pertaining  to  continued  fractions.  Section  A. 2  introduces  and  proves  several  results  about 
binary  quadratic  forms  and  describes  their  connection  to  continued  fractions.  Section  A. 3 
analyzes  how  ideals  are  related  to  quadratic  forms  and  section  A. 4  analyzes  how  lattices  are 
closely  tied  in  with  ideals.  Section  A. 5  uses  the  connection  between  quadratic  forms  and 
lattices  to  prove  the  general  distance  formula,  which  is  used  in  section  A. 6  to  explain  why 
SQUFOF  works. 

A.l  Continued  Fractions 

One  tool  used  by  many  different  algorithms  is  the  continued  fraction  expression  for  \/N , 
where  N  is  the  number  to  be  factored.  This  expression  is  calculated  recursively [10]: 


x0  =  Vn,  bo  =  L^oJ 

(A.l) 

7  i  h  =  bd 

(A.2) 

Xi—i  b{— i 

Vn  =  b0  +  , 

(A.3) 

'i  + 


b2+... 
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Observe  that  solving  equation  (A. 2)  for  Xi-i  gives  =  A-i  +  Repeatedly  substi¬ 
tuting  this  into  itself  gives  equation  (A. 3). 

Before  developing  the  theory  too  much  further,  allow  me  to  offer  one  example  of  how 
this  works,  just  so  that  the  pieces  make  sense  to  you.  To  simplify  the  expansion  some,  the 
integers  taken  out  in  the  third  step  of  each  line  are  the  bp. 


Example  9. 


Xi  = 


x2  = 


%3  = 


x0  =  Vil,  bo  =  6 

1  _  yiI+6  _  o  I  yTl— 4 

yiI-6  5  _l"  5 

5  _  QiI+4  _  9  I  vTL-6 

vTI-4  5  5 

5  _  VIl+6  _  19  1  x/TT— 6 

vTT-6  “  1  ~~  iZ  "l"  1 


^  =  6+  A 


v/4l  =  6  + 


2+  — 
^2 


v/4l  =  6  + 

"s/41  =  6  +  - 


2+ 


2+  — 


2+- 


12+^ 


From  this  step,  it  is  evident  that  aq  =  x\  and  the  cycle  repeats  from  here.  Observe  that 
if  the  sequence  starts  with  x0  =  \/41  +  6,  then  also  x3  =  x0-  Regardless,  the  expression  on 
the  right,  if  truncated  at  any  point  provides  a  rational  approximation  to  ^41.  Often  this 
will  be  written  as  merely  [6,  2,  2, 12, ...]  to  save  space  or  [6,  2,  2, 12]  to  indicate  that  this  part 
repeats.  The  various  numbers  on  the  left  have  some  important  properties  that  we  will  now 
analyze  in  some  depth. 

Throughout,  assume  that  N  is  an  odd  positive  integer  and  is  not  a  perfect  square.  For 
number  theory  purposes,  let 


x0  = 


Vn  +  -P-i 
Qo 


where  P__i,Q0  are  integers  chosen  such  that 


P\  =  N  (mod  Qo),  0  <  P_i  <  y/N,  and  \VN  -  Q0\  <  P^.  (A.4) 


There  are  many  ways  of  doing  this1.  The  recursive  formulas  are: 

1Choosing  Xq  =  \/N  +  f\/iVJ ,  so  that  P_i  =  fx/iV J  and  Q0  =  1  is  one  possibility.  Choosing  Xq  = 
where  P-i  =  f\/AtJ  or  —  1,  such  that  P_i  is  odd,  is  another  possibility. 
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*£*+ i 


Xi  -  bi 


h  =  L^J ,  i  >  0 


Formally,  the  assumed  equation  is: 


•^i+1 


Qi 


VN  +  Pi  ,  ,  VN-Pi+1  ,^n 

=  bi+ 1  H - — - ,  %  >  0 


(A.5) 


VF  -P»  Qi+i  “lTi  '  Qt+i 

Note  that  this  equation  serves  as  a  definition  of  Qi,  Pt,  Qi+ 1,  Pi+1  e  Q,  so  that  these  equations 
are  true  regardless  of  the  conditions  on  these  variables.  Theorem  A. 1.1  provides  some  well- 
known  fundamental  properties  and  identities  of  continued  fractions.  In  [18],  Hans  Riesel 
provides  very  clear  proofs  of  most  of  this. 

Theorem  A. 1.1.  [18]  In  the  continued  fraction  expansion  of  xq  satisfying  (A. 4),  each  Xi 
reduces  to  the  form  ^ P'~1 ,  with  (a)  N  =  Pf  +  QiQi+i,  (b)  Pi  =  biQi  —  Pi- 1,  (c)  bi  = 
Pl~ 1  —  (d)  0  <  Pi  <  VN,  (e)  \VN  —  QA  <  Pi_ i,  (f)  Qi  is  an  integer,  and  (g) 

Qi+i  =  Qi- 1  +  bi(Pt- 1  —  Pi).  Furthermore,  (h)  this  sequence  is  eventually  periodic. 


Proof: 


(a)  From  (A.5),  the  equation  requires  that  N  —  Pf  +  QiQi+1. 

(b)  It  is  evident  from  simplifying  the  expression  on  the  far  right  of  (A.5)  that 


a/F  +  Pi  \/F  +  bi+\Qi+ 1  —  Pi 


i+ 1 


Qi+ 1  Qi+ 1 

Therefore,  Pi+1  =  bi+iQi+1  -  Pt. 

(c)  For  i  =  0,  by  the  assumption  |  \/F  —  Qo\  <  P- 1 

Qo  <  VN  +  P- 1 


Therefore, 


5n  = 


V^  +  P-l 

Qo 


>  l 


For  i  >  0,  bi-i  =  \_Xi- ij.  By  the  definition  of  floor,  Xi-i  —  1  <  bi- 1  <  Xi-\.  If  bi- 1  =  Xj_i, 
then  the  continued  fraction  [60, 6]_,  is  rational  and  is  equal  to  x0  =  ,  which 
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is  irrational  since  N  is  not  a  perfect  square.  Therefore,  Xj_i  —  1  <  i  <  so  that 

0  <  Xi_ i  —  &j_i  <  1.  Therefore,  Xi  =  i  >  1,  so  that  bi  =  >  1. 

Note  that  there  is  no  integer  between  [v^/Vj  +  Pi- 1  and  \fN  +  Pi- 1,  so  it  is  trivial  that 


pN+Pi- 1 

VA'  4/’.  i 

Qi 

Qi 

(d-e)  The  statements  \VN  —  Q{\  <  P,_i  and  0  <  P,_i  <  \fN  may  be  proven  inductively. 

Base  case:  i  =  1 

P0  =  [VN\  or  [\/N\  —  1,  so  by  definition  0  <  P0  <  \/N. 

Since  x0  meets  (A. 4),  j \/N  —  Q0 \  <  P_i. 

Induction:  Assume  \y/N  —  Qi\  <  P*_i  and  0  <  P,_i  <  \/N . 

Note  that  these  assumptions  require  that  0  <  Qi  <  2\/iV.  From  (c),  0  <  Xi  —  bi  <  1 

means  0  <  <  1.  Since  Qi  >  0,  0  <  \/N  —  Pi  <  Qi.  From  the  left  side  of  this, 

Pi  <  V ~N .  Now,  either  Qi  <  \J ~N  or  Q{  >  a /N. 

Case  1:  If  Qi  <  VN,  then  \/N  —  Pi  <  Qi  <  VN,  so  that  p  >  0. 

Case  2:  If  Qi  >  \/N,  then  by  (b),  Pi  =  biQi  —  P*_i  >  bi\/N  —  \/N  =  ( bi  —  1)a /N  >  0. 

Therefore,  0  <  P*  <  \zrN. 

Since  xi+\  >  1,  it  is  trivial  that  Qi+ 1  <  \fN  +  P,:  so  that  showing  |  y/N  —  Qi+ 1|  <  Pj 
reduces  to  showing  Qi+ \  >  \fN  —  Pi.  Since  1  =  ^ nPi  =  this  is  equivalent  to 

showing: 


VN  +  Pt 
Qi 


>  l. 


(A.6) 


Assume  the  contrary,  that  Qi  >  \QN  +  P*.  Then, 


bi(VN  +  Pi)  -  Pi  <  biQi  -  Pi  =  Pi _i  <  VN, 

biVN  +  Pi(bi  -  1)  <  VNt 

VN(bi  -  1)  +  Piibi  -  1)  <  0, 

{h  -  l)(VN  +  Pi)  <  0. 

But  VN  and  Pi  are  positive,  so  this  implies  bi  <  1,  contradicting  Theorem  A. 1.1  (c).  There- 
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fore,  (A. 6)  holds. 

(f)  The  fact  that  N  —  Pf  +  QiQi+i  requires  that  Qi+ 1  =  N  QPi  .  In  order  to  show  that 
V*  Qi  is  an  integer,  the  statements  that  Qt  is  an  integer  and  Qi  \  N  —  Pf  may  be  proven 
inductively. 

Base  case:  i  =  0 

By  definition,  Q0  is  an  integer  and  Pf1  =  N  (mod  Q0).  But  P0  =  P-i  (mod  Q0 ),  so 
P02  =  N  (mod  Q o).  Therefore,  Q0  \  N  —  P02. 

Induction:  Assume  for  some  i,  Qi  is  an  integer  and  Qi  \  (N  —  Pf).  Then,  since  N  = 
Pf  +  QiQi+i,  Qi+ i  =  so  that  since  Qi  \  (N  -  Pf),  Qi+1  is  an  integer.  Also,  Qt  = 

Nq.+)  ,  so  that  since  Qi  is  an  integer,  Qi+\  \  (N  —  Pf),  so  that  Pf  =  N  (mod  Qi+ 1).  Since 
Pi+ i  =  bi+1Qi+1  -  Pi,  Pi+i  =  Pi  (mod  Qi+i).  Therefore,  Pf+1  =  N  (mod  Qi+i),  so  that 
Qi+i  |  (N  —  Pf+1)  and  the  induction  is  complete. 

(g)  Solving  (b)  for  bi  gives  Pi~T+Pi  =  6j.  Multiply  by  (P*_ i  —  Pf)  to  obtain: 


Pf-i-Pf 

Qi 


Rearranging  and  adding  ff  gives: 


bi(Pi~  i  —  Pf 


N  -Pf 
Qi 


N 


p  2 

ri- 1 


Qi 


+  bi(Pi- 1  —  Pf 


Qi+l  —  Qi- 1  +  bi(Pi- 1  —  Pi) 


(h)  Since  each  a:*  and  thus  the  entire  sequence  that  follows  it  is  defined  by  the  two  integers 
Qi  and  P,_i,  limited  by  the  bounds  0  <  Qi  <  2 y/N  and  0  <  P*  <  VN,  there  is  only  a  finite 
number  of  distinct  xf  s.  Therefore,  for  some  tt  and  some  k,  Vi  >  k  Xi  =  Xi+7T.  QED 

The  fact  that  each  Xi  reduces  to  the  form  is  important  for  computational  effi¬ 

ciency  because  this  together  with  (c)  imply  that  floating  point  arithmetic  is  not  necessary  for 
any  of  these  calculations.  Also,  by  use  of  (b)  and  (g),  the  arithmetic  used  in  this  recursion 
is  on  integers  <  2  y/N. 
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One  application  of  continued  fractions  is  rational  approximations. 


VN  =  b0  + 


1 


&i  + 


i 

62+... 


If  this  continued  fraction  is  truncated  at  any  point,  the  result  is  an  approximation  to 
y/N.  One  might  imagine  that  it  is  necessary  to  start  simplifying  at  the  lower  right  end  of 
this  expression  to  obtain  this  approximation.  However,  Theorem  A.  1.2,  also  included  in  [18], 
provides  a  simpler  answer. 


Theorem  A.  1.2.  Let: 


A- 1  —  1,  A0  —  bo,  Ai  —  biAi_i  +  Aj_2,  i  >  0 
B- 1  =  0,  B0  =  1,  Bi  =  bjBi_  1  +  Bi_ 2,  i  >  0 


Then  for  i  >  0,  [60,  0,  ■■■bi]  =  jf  arid  A2_x  —  B2_  1 N  =  (— 1  )lQi. 

Note  that  the  last  equation  gives  A2_x  =  (— 1  )lQj  (mod  N ).  Although  this  equation 
will  change  some  when  generalized  to  other  continued  fractions,  these  denominators  {Qi} 
will  consistently  be  referred  to  as  pseudo-squares.  A  proof  of  this  theorem  is  given  in  [18]. 
Therefore,  instead  of  reproducing  the  proof,  I  will  provide  an  example: 

Example  10. 

x0  =  \/403  +  20  =  40  +  \At03  -  20 

_  1  _  x/403+20  _  10  I  VT03-19 

''  1  “  ^403-20  “  3  -  10  "T  3 

_  3  _  n/403+19  _  o  1  V403-9 

2  VT03-19  14  _l"  14 

_  14  VT03+9  _  1  1  yTol-M 

3  VT03-9  23  1  23 

„  _  23  _  n/403+14  _  q  j.  7403-13 

X4  “  x/403-14  “  9  “  °  ~r  9 

9  y/403+13  1  ,  7403-13 

5  7403-13  26  26 

_  26  _  7403+13  _  q  ,  7403-14 

X6  “  7403-13  “  9  -  0  9 

_  9  _  7403+14  _  i  ,  7403-9 

7  7403—14  23  ”1'  23 

From  this,  a  table  may  be  used  to  recursively  calculate  the  approximation2: 


2Actually,  in  this  case  b0  =  40,  but  since  that  would  be  an  approximation  to  x0  =  V403+-  20,  subtracting 
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i  -1  0  1  2  3  4  5  6 

bi  20  13  2  1  3  1  3 

Ai  1  20  261  542  803  2951  3754  14213 

Bi  0  1  13  27  40  147  187  708 

filling  in  and  B \  from  left  to  right.  From  the  last  column,  ^403  ~  14213/708. 

Since  the  continued  fraction  is  eventually  periodic,  it  is  reasonable  to  consider  that  when  it 
loops  around  on  itself,  the  terms  being  considered  may  have  come  from  some  terms  “earlier” 
in  the  recursion.  Example  10  provides  some  indication  as  to  how  the  recursive  formulas 
may  be  reversed,  as  {Qi\  and  {bi}  are  symmetric  about  x$,  so  that  after  x$  these  numbers 
are  cycled  through  in  reverse  order.  Lemma  A.  1.3  addresses  how  each  bi  is  calculated  two 
different  ways  and  Lemma  A.  1.4  shows  that  by  exchanging  these  two  related  expressions, 
the  direction  is  reversed. 

Lemma  A. 1.3. 

V~N  +  Pi  .  ,  y/N  +  Pi-i  . 

1  = 


Proof:  The  second  part  of  this  equation,  that  \_^QPi  =  bt  follows  from  the  definition 
of  bi. 

Theorem  A. 1.1  (e)  implies  that  Qi  >  y/N  —  P,_i.  Therefore, 


VN  +  Pt 
Qi 


\ f~N  +  biQi  —  Pi-i 

Qi 


bi  + 


Vn  -  p^ i 

Qi 


=  bi.  QED 


Considering  Example  10,  it  is  then  natural  to  suspect  that  the  mechanism  for  going  in 
the  opposite  direction  will  be  precisely  the  same  as  the  standard  approach,  except  that  the 
numerator  is  changed  first.  Note  that  this  same  change  (with  the  exception  of  Co)  could  be 
achieved  by  merely  changing  the  sign  of  Pt- 1 . 


20  from  bo  yields  an  approximation  to  \/403. 
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Lemma  A. 1.4.  Let  x^b^P^Qi,  and  N  be  as  in  Theorem  A. 1.1,  %  >  0.  Let  yo 
and  let  cq  =  |j/oJ  ■  Define  inductively  yj  =  —  - .  Then  cq  =  bl+\  and  yj  = 


_=  VN+Pj+i 

Qi-\- 1 

y/~N-\-Pj—j+ 1 
Qz — j+1  ^ 


J  >  0. 


Proof:  By  (A. 6)  and  Lemma  1,  Co  =  |_2/oJ  =  |_  V^q+^i+1  J  =  ^?:+ 1  •  By  mathematical 
induction  it  suffices  to  prove  the  case  j  =  1.  Using  Theorem  A. 1.1 


1 

Ih  = - 

Vo  ~  c0 


VN+Pi+1 
Qi-\- 1 


-h 


'i+ 1 


1 

yA +Pj+i— bj+iQi+i 
Qi+l 


1  =  v^  +  P, 

VN-Pi  N~pj 

Qi+l  Qi+l 

This  demonstrates  an  important  fact  about  continued  fractions,  the  fact  that  the  direction 
of  the  sequences  of  pseudo-squares  and  residues  can  be  reversed  (i.e.  the  indices  decrease) 
by  making  a  slight  change  and  applying  the  same  recursive  mechanism. 

Using  Lemma  A. 1.4,  x3  may  be  used,  for  example,  to  find  X2  and  x\.  Continuing  this 
process,  denote  the  terms  before  Xo  as  X-±,X-2, ....  Define  Q-i  and  P_j  similarly3.  Example 
11  demonstrates  this  with  the  continued  fractions  from  Example  10: 

Example  11.  x3  =  v/493+9  and  P3  =  14,  so  let  y0  =  v/^3+14  to  obtain 


y/N  +  Pj 
Qi 


QED 


x/403+14  1  , 

y°  23  1  23 

_  23  _  vToS+Q  _  0  1  y/403-19 

yi  VT03-9  14 


14 


_  14  _  VT03+19  _  IQ  I  VT03 

V2  -  ^403-19  “  3  -  +  3 


-20 


y3  = 


y/403-20 


_  y/403+20  _  _|_  t/403 


VT03-20 


_  1  _  x/403+20  _  iq  1  VT03-19 

VA  ~  VU03-20  ~  3  “  40  +  3 


7,  =  3  VT03+19  _  9  1 

y5  VT03-19  14  "r  14 


Then,  just  as  y2  gives  xx  =  x/49|+2°,  y4  gives  x_i  =  %/49|+19  and  y5  gives  x_2  = 
Combining  periodicity  with  reversibility  strengthens  Theorem  A.  1.1  (h). 


3Since  j/0  meets  (A. 4)  and  the  same  recursive  formula  is  applied,  it  is  clear  that  Theorem  A. 1.1  still 
applies  to  negative  indices. 
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Lemma  A. 1.5.  There  exists  a  positive  integer  n  such  that  Vi  x*  =  xi+7r7  i  not  necessarily 
positive. 

Proof:  From  the  proof  of  Theorem  A.  1.1  (h)  there  are  k  and  n  such  that  Vi  >  kj 
Xi  =  xi+7T.  Essentially,  this  is  equivalent  to  proving  that  there  is  no  lower  bound  for  k. 
Assume  the  contrary,  that  there  is  some  lower  bound  k.  Let  k  and  n  be  the  smallest  such 
integers.  Then  x k  =  Xfc+7r.  But  by  Lemma  A.  1.4  Xk-i  =  Xk+n-i,  so  that  k  —  1  also  meets  this 
criteria,  violating  the  assumption  that  k  is  the  smallest  such  integer.  Therefore,  Vi  Xi  =  xi+1T. 
QED 

Throughout,  n  will  consistently  denote  the  period,  even  when  considering  this  period  in 
the  context  of  quadratic  forms  or  lattices. 

Often  the  continued  fraction  may  have  other  characteristics  that  are  interesting  besides 
its  periodicity.  For  factorization,  continued  fractions  with  symmetries,  such  as  at  Xo  and  X5 
from  Example  10,  will  be  especially  important.  If  the  starting  condition  near  some  point  is 
the  same  in  both  directions,  the  entire  sequence  will  be  symmetric  about  that  point.  This  is 
the  point  of  Lemma  A.  1.6. 

Lemma  A. 1.6.  Let  x0  =  meet  (A. 4)  such  that  Qo  \  2P_i.  The  sequence  of  pseudo¬ 

squares  is  symmetric  about  Qo,  so  that  \/i  Q,  =  Q-i. 

Proof: 

Observe  that  0  <  VN  —  P0  <  Q0,  with  P0  =  b0Q0  —  P_i ,  so  that 


0  <  \/~N  —  b0Q0  +  P—i  <  Qo ■ 


There  can  only  be  one  possible  integer  value  of  b0  that  satisfies  this  inequality.  Since  0  < 


Vn  -  P_! 
Let  y- 1 


<  Qo,  b0  =  2P_i/Qo  satisfies  this  inequality,  so  that  P0  =  P_i. 
=  ^q/1-  Then>  bV  Lemma  2>  Vo  =  =  =  xo 


Therefore,  the  sequence  of  pseudo-squares  will  be  symmetric  about  Qo,  since  in  either 


direction  the  first  continued  fraction  term  is  the  same.  Therefore,  Qi  =  Q-i-  QED 
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The  presence  of  one  point  of  symmetry  allows  a  proof  that  another  point  of  symmetry 
exists  and  that  a  factorization  of  N  may  be  obtained  from  this  symmetry4: 

Theorem  A. 1.7.  Let  s  =  |_§  J .  where  tt  is  the  period  from  Lemma  A. 1.5.  If  tt  is  even, 
Vi  Qs+i  =  Qs-i,  but  Qs  7^  Q0  and  Qs  |  2 N.  If  n  is  odd,  Vi  <3s-h+i  =  Qs-i  and  either 
gcd (QS,N)  is  a  nontrivial  factor  of  N  or  —1  is  a  quadratic  residue  of  N. 

Proof: 

Case  1:  If  tt  is  even,  n  =  2s.  Then,  by  Lemmas  A. 1.6  and  A. 1.5,  Qs+i  =  Q-s-i  = 
Q2s-s-i  =  Qs-i •  Since  Qs+i  =  and  Qs_  1  =  -  this  simplifies  to  Pf  =  Pf_i,  but 
since  Vi  Pi  >  0,  this  provides  Ps  =  Ps-\. 

Now  Qs  =  P3+^3~1  =  so  that  Qs  |  2 Ps. 

Assume  Qs  =  Qq.  If  Qs  is  even,  then  P0  =  Ps  =  1  (mod  2)  and  if  Qs  is  odd,  Qs  \  Ps. 
Either  way,  there  is  then  a  unique  integer  in  the  range  (y/N  —  Q o,  y/N)  satisfying  these 
conditions,  so  that  Ps  =  Pq.  Therefore,  xs  =  =  Xo,  contradicting  the  fact  that  tt  is 

the  smallest  positive  integer  such  that  Vi  Qi  =  Qi+n.  Therefore,  Qs  ^  Qo- 

Now  N  =  Pf  +  QsQs+i,  so  it  is  apparent  that  if  Qs  is  odd,  then  Qs  \  Ps,  so  that  Qs  \  N. 
Conversely,  if  Qs  is  even,  then  (Qs/ 2)  |  Ps,  so  that  (Qs/ 2)  j  N.  Either  way,  Qs  \  2 N . 

Case  2:  If  tt  is  odd,  tt  —  2s  +  1.  Then,  by  Lemma  A. 1.6  and  A. 1.5,  Qs+i+ 1  —  Q-s-i- i  = 
Q2s+l— s— i—  1  Qs—i- 

Specifically,  Qs  =  Qs+i,  so  that  N  =  Pf  +  QSQS+1  =  Pf  +  Q2S,  so  that  Pf+l  =  —  Q2S 
(mod  N ).  If  gcd(Qs,  N)  >  1,  this  is  a  nontrivial  factor  of  N,  and  the  proof  is  done.  Therefore, 
assume  that  Qs  and  N  are  relatively  prime,  so  that  Qj1  (mod  N )  exists.  Then  (QL1)2^^!  = 
— 1  (mod  N ).  Then  Q^Pg+i  is  a  square  root  of  —1  modulo  N.  QED 

One  final  concept  that  will  appear  much  more  important  in  later  sections  is  equivalence. 
Define  the  set  T  to  be  set  of  all  numbers  of  the  form  ^+F  such  that: 

P2  =  N  (mod  Q),  (A. 7) 

4This  was  actually  discovered  in  the  opposite  order.  It  was  clear  that  ambiguous  forms  that  met  this 
criterion  provided  a  factorization  but  was  later  realized  that  these  same  forms  produced  symmetry  points. 
This  was  first  noticed  by  Gauss  [8]  and  first  applied  by  Shanks  [23]. 
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Then  define 


r  =  {ieT:0<P<  Vn,  I  Vn-q  I  <  P}. 

An  element  x  G  T  is  reduced  if  x  G  T*.  For  x,?/  6  F,  i  is  equivalent  to  y  if  x  appears 
in  the  same  continued  fraction  expansion  as  y  and  it  is  trivial  that  this  is  an  equivalence 
relation  on  T*.  Extending  this  to  all  of  T  requires  a  lemma  relating  elements  of  T  —  T*  with 
elements  of  T*. 

Lemma  A. 1.8.  Let  x  =  x0  e  T  —  T*.  x0  may  be  reduced  by  a-p-plying 


xi+ 1  =  —  .  ,  h  =  [xi  -  1/2J,  i  >  0  (A. 8) 

Xi  —  bi 

until  \Qi\  <  2 y/N  for  some  i  and  then  applying  equation  (A. 2)  normally  until  Xk  G  T*  for 
some  k  >  0. 


Proof:  The  choice  of  b,  yields  that 
\\Qi\  +  \/N.  Therefore, 


VN-Pj 

Qi 


<  |  after  the  first  step,  so  that  \Pt\  < 


1 

N-Pf 

1  ~ 

Qi 

sfN-Pi 


Qi  ||V^+^I 

<  |(2 y/N+±\Qi\ 

=  VN  +  \\Qi\ 


so  that  \Qi\  will  decrease  as  long  as  \Qf  >  | y/N.  When  \Qr\  <  2 y/N,  revert  back  to  the 
standard  formula  for  br.  There  are  three  cases  for  what  Qr  is: 

Case  1:  0  <  Qr  <  y/N .  In  this  case,  it  is  clear  that  xr+\  will  be  reduced. 

Case  2:  y/N  <  Qr  <  2 y/N.  In  this  case,  if  Pr  >  0,  then  xr+\  will  be  reduced.  Otherwise, 
|  Pr  |  <  y/N,  so  that  xr+\  will  be  in  Case  1. 

Case  3:  —2 y/N  <  Qr  <  0.  Then  y/N  <  Pr  <  y/N  +  \Qr\,  yielding  Qr+\  <  2 y/N  +  \Qr\. 
If  Qr+i  <  2 y/N,  it  is  in  Case  1  or  Case  2.  If  2 y/N  <  Qr+\  <  2 y/N  +  \Qr\,  the  choice  of  br 
provides  y/N  +  Pr  >  Qr+i,  so  that  0  <  Pr+\  <  y/N,  so  that  xr+2  will  be  in  Case  1.  QED 
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I  will  provide  one  example  of  reduction: 


Example  12. 

„  _  y/403+267  _  _9  I  \/403-l 
-134  Z  ”r  -134 

„  _  -134  _  V403+1  _  o  |  v/403-23 

Xl  ~  \/403— 1  “  -3  “  °  -3 

-3  y/403+23  i  ,  y/403— 19 

X2  v/403-23  42  “l"  42 

_  42  _  \/403+19  _  Qr)  .  VT03-20 

X3  “  \/4C)3— 19  “  1  ~~  1 

Lemma  A.  1.8  defines  a  map  from  T  to  T*  (Elements  of  T*  are  mapped  to  themselves). 
Then  two  elements  are  equivalent  if  their  corresponding  elements  of  T*  are  equivalent  and  it 
is  clear  that  this  is  still  an  equivalence  relation.  Essentially,  this  equates  to  saying  that  two 
numbers  x  and  y  are  equivalent  if  their  continued  fraction  expansions  have  the  same  “tail” , 
so  that  after  a  certain  number  of  terms  of  each  they  have  identical  cycles.  Appendix  A. 2 
will  define  a  different  equivalence  relation  on  binary  quadratic  forms  and  then  prove  that  it 
corresponds  to  this. 


A. 2  Binary  Quadratic  Forms 


A  fuller  account  of  binary  quadratic  forms  can  be  found  in  Gauss’s  [8]  and  in  Buell’s  [2], 
However,  here  are  the  necessary  fundamental  ideas. 

A  binary  quadratic  form  is  a  polynomial  of  the  form  F(x,  y )  =  ax2-\-bxy+cy2  ,x,  y,  a,  b,c  G 
Z  (Often  this  is  abbreviated  as  ( a,b,c )).  In  some  sense  then,  a  quadratic  form  may  be 
considered  to  be  the  set  of  all  the  numbers  it  can  represent  for  various  values  of  x  and  y. 
Thus,  two  quadratic  forms  are  equivalent  if  they  represent  the  same  set  of  integers.  It  is 
evident  that  if  one  form  is  transformed  into  another  by  the  substitution 


X 

a  b 

x' 

y 

c  d 

y' 

ad  —  be 


±1, 


(A.9) 


then,  since  this  matrix  is  invertible,  the  two  forms  are  equivalent.  As  one  further  useful 
distinction,  (A.9)  is  proper  if  its  determinant  is  +1  and  improper  if  its  determinant  is  —  1. 
The  symbol  (~)  will  only  apply  to  proper  equivalence.  For  the  purpose  of  factorization,  the 
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interesting  forms  are  those  that  can  be  improperly  transformed  into  themselves,  referred  to 
as  ambiguous  forms. 

Example  13.  F(x,y )  =  —14a;2  +  10 xy  +  5 y2  is  transformed  into  itself  by  the  substitution: 


so  (—14, 10,  5)  is  ambiguous. 

Denote  the  set  of  all  quadratic  forms  with  discriminant  A  by  F^,  or  often  just  F.  The 
next  obvious  question  is  the  organization  of  all  of  the  quadratic  forms  equivalent  to  some 
given  form.  Since  there  are  an  infinite  number  of  forms  equivalent  to  any  form,  the  search 
must  be  narrowed  some  by  first  defining  reduced  forms. 

Definition  A. 2.1.  A  quadratic  form  ax2  +  bxy  +  cy2,  with  positive  discriminant  A  =  62  —  4ac 
is  reduced  if: 

0  <b<VA  (A. 10) 

|\/A  —  2|a||  <  b  (A. 11) 

Note  that  A  =  b2  —  4 ac  and  (A.  10)  require  that  ac  <  0,  so  that  a  and  c  must  have 
opposite  signs. 

Making  Gauss’s  description  of  the  organization  make  some  sense  will  require  one  more 
of  his  definitions: 

Definition  A. 2. 2.  Two  forms  F(x,y )  =  ax 2  +  bxy  +  cy 2  and  F'(x,y )  =  a'x 2  +  b’xy  +  c'y2 
are  adjacent  if  c  =  a'. 

To  each  quadratic  form,  there  is  a  unique  reduced  equivalent  form  adjacent  to  it  on 

each  side5,  and  since  (A.10-A.11)  imply  a  finite  number  of  possible  coefficients,  this  process 

5 Although  Gauss  had  a  recursive  mechanism  for  finding  these,  continued  fractions  provide  a  sufficient 
mechanism  for  this  that  will  be  defined  momentarily.  Note  that  reversal  suddenly  becomes  trivial. 
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eventually  repeats,  forming  a  cycle.  The  important  aspect  of  this  is  that  the  cycle  is  actually 
all  of  the  reduced  forms  equivalent  to  the  first  form: 

Theorem  A. 2. 3.  [8]  If  the  reduced  forms  F,F'  are  properly  equivalent,  each  of  them  will  be 
contained  in  the  period  of  the  other. 

Gauss  proves  this  in  Article  193  of  [8],  Lenstra  proves  this  in  [12],  and  it  is  a  corollary  of 
Lemma  A. 4. 7  in  Appendix  A. 4.  Therefore  the  proof  is  omitted. 

For  now,  the  important  detail  is  that  the  quadratic  forms  correspond  directly  to  the 
elements  of  T,  and  that  the  reduced  quadratic  forms  correspond  to  elements  of  T*.  Note 
that  the  elements  of  T  have  attached  indices,  where  the  important  trait  of  the  index  is 
whether  it  is  odd  or  even.  Define  a  map  from  T  to  F  by 


$t,f  :  T  — >  F 


Vn-p \ 
Qi 


F,(x,y)  =  Qi(- l)‘x2  +  2P,xy  +  Qi+i(-l)‘+1s2 


(A. 12) 


The  inverse  map  is 


<f>F,T  :  F  — >  T 

2  ;  2  ™  (A-12  ) 
ax 2  +  bxy  +  cy2  — >  ay  =  ^ - G  T 

where  the  discriminant  of  the  quadratic  form  is  A  is  the  element  of  T  is  given  either  an 
even  or  odd  index  as  a  is  positive  or  negative,  respectively.  Note  that  A  =  b2  —  4 ac  gives 
4a  |  A  -  b2,  so  that  Xi  really  is  in  T. 

§A.l  defined  an  equivalence  on  T  and  suggested  that  it  corresponded  with  an  equivalence 
of  binary  quadratic  forms.  Theorem  A. 2. 4  formalizes  this: 

Theorem  A. 2. 4.  Under  the  mapping  d^F.,  the  equivalence  classes  of  T  correspond  to  the 
equivalence  classes  of  F.  That  is,  for  Xi,Xj  G  T  corresponding  to  Fi  =  Tn-y  (ay),  iy  = 
$T,F(^j)  G  F,  respectively,  ay  ~  ay  if  and  only  if  Fi  ~  Fj. 

Proof: 

Let  Xi  ~  Xj.  Since  ay  and  ay  must  be  in  the  same  continued  fraction  expansion,  assume 
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without  loss  of  generality  that  j  —  i  +  1.  The  other  cases  may  be  easily  derived  from  this 
case.  Then  the  quadratic  form  related  to  xt  is  given  in  (A.  12).  The  substitution 


X 

0 

l 

x' 

y 

-l 

(-i)Vi 

y' 

transforms  Ft  into 

Qi+i{—iy+1x2  +  2  Pi+ixy  +  Qi+2(—l)l+~y2 

Observe  that  this  matrix  has  determinant  1,  so  that  this  equivalence  is  proper. 

In  order  to  prove  the  converse,  that  the  x^s  related  to  equivalent  quadratic  forms  are 
equivalent,  by  Theorem  A. 2. 3,  observe  that  the  last  coefficient  of  the  quadratic  form  related  to 
Xi,  Qi+ 1(— 1)*+1,  is  then  the  first  coefficient  of  the  quadratic  form  related  to  xl+\.  Therefore, 
these  two  forms  are  adjacent  and  thus  equivalent.  QED 

The  real  value  of  quadratic  forms  is  the  composition  of  quadratic  forms.  In  Article  236 
of  [8],  Gauss  provides  a  very  flexible  definition  of  composition.  Gauss  defines  composition 
as  multiplying  two  quadratic  forms  together  and  then  making  a  substitution  to  simplify  this 
into  another  binary  quadratic  form.  The  algorithm  he  provides  is  very  complicated,  allowing 
for  choices  of  variables  along  the  way  that  permit  the  result  to  be  any  quadratic  form  in  the 
resulting  equivalence  class.  The  result  of  composition  should  be  predictable,  so  definition 
needs  to  be  limited  some.  Shanks  and  Buell  both  provide  a  significant  simplification  of  this 
algorithm.  The  symbol  (*)  will  consistently  be  used  for  composition. 

Proposition  2.  [2]  Let  Fx  =  (op,  bi,ci)  and  F2  =  (a2,  b2,c2 )  be  primitive  forms  of  discriminants 
d\  and  d2,  respectively,  such  that  d\  =  A n\  and  d2  =  A n\  for  integers  ri\  and  n2  and  A, 
with  A  =  gcd(di,  d2).  Let 


m  =  gcd(airi2,  a2ni, 


bin2  +  b2ni 
2 


Then  the  congruences 
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mriiB  =  mbi  (mod  2ai) 
rrir^B  =  mb2  (mod  2a2) 

rn(b\Ti2  +  b2ni)B  =  771(6162  +  Anin2)  (mod  4aia2) 
are  simultaneously  solvable6  for  an  integer  B ,  and  the  composition  of  Fl  and  F2  is: 


of  discriminant  A. 


Fi*  F2  = 


,d\a2  ( B 2 

,  2  ’ 
rnz 


A  )m2 


4oi(22 


See  [21]  for  a  derivation  of  this  in  the  case  where  the  discriminants  are  equal  or  [2]  for  a 
proof  of  this  case.  Buell  [2]  also  provides  the  substitutions  that  would  be  needed  for  Gauss’s 
definition  of  composition. 

The  next  question  is  how  this  operation  is  related  to  equivalence. 


Theorem  A. 2. 5.  If  F\  ~  F2,  then  F  *  F\  ~  F  *  F2. 

Gauss  proves  this  in  Article  237-239  of  [8]. 

Therefore,  composition  treats  the  equivalence  classes  in  the  convenient  manner.  These 
equivalence  classes  are  then  the  elements  of  the  class  group,  with  composition  as  the  group 
operation.  The  application  of  Theorem  A. 2. 5  is  that  it  doesn’t  matter  which  form  is  used  to 
represent  an  equivalence  class. 

The  significance  of  ambiguous  forms  for  factorization  has  been  mentioned  some  above.  It 
is  evident  that  if  one  form  is  ambiguous,  then  its  entire  equivalence  class  is  also  ambiguous. 
Lemmas  A. 2. 6  generalizes  the  reasons  to  be  interested  in  these  classes. 


Lemma  A. 2. 6.  An  ambiguous  equivalence  class  contains  two  points  of  symmetry,  that  is, 
pairs  of  reduced  adjacent  forms,  (c,  6,  a )  and  (a,  6,  c)  in  the  cycle  that  are  the  symmetric 
reverse  of  each  other.  Let  a  be  the  connecting  term  of  either  symmetry  point.  Either  a 
divides  the  determinant,  ora/2  divides  the  determinant. 

Proof: 


6By  convention,  choose  the  answer  with  the  smallest  absolute  value. 
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Let  A  be  an  ambiguous  equivalence  class  and  let  F  =  ax 2  +  bxy  +  cy2  G  A.  Let  F'  = 
ex2  +  bxy  +  ay2.  Then  since  F  G  A,  there  is  a  substitution  of  determinant  —1  that  maps  F 
into  itself.  Since  the  obvious  substitution  to  exchange  x  and  y  in  F  has  determinant  —1,  the 
product  of  these  two  is  a  proper  substitution  that  transforms  F  into  F' .  Therefore,  F'  G  A, 
so  that  if  F  is  Fq,  F'  is  Fj  for  some  j.  Then  that  F\  must  be  the  reverse  of  Fj_ i,  and 
so  forth.  Now,  if  j  is  even,  then  by  this  process  Fj/2  is  its  own  reverse.  However,  by  the 
definition  of  being  reduced,  the  end  coefficients  of  each  form  must  have  opposite  sign,  so  this 

is  impossible.  Therefore,  j  must  be  odd,  and  then  Fj-i  is  the  reverse  of  Fj+i . 

2  2 

At  this  point,  observe  that  since  the  end-coefficients  alternate  signs,  the  entire  period 
must  be  even.  By  the  same  arguments  as  Theorem  A. 1.7,  one  could  show  that  there  must 
be  another  point  of  symmetry  with  the  property  that  Vi  Qs+i  =  Qs-i ,  but  such  that  Qs  is 
not  the  same  as  the  connecting  term  at  the  first  symmetry  point.  The  two  quadratic  forms 
containing  Qs  as  an  end  coefficient  then  meet  the  criteria. 

The  fact  that  either  a  divides  the  determinant,  or  a/2  divides  the  determinant  was  proven 
in  Theorem  A. 1.7,  since  the  determinant  is  41V.  QED. 

Note  that  Theorem  A. 1.7  described  two  different  types  of  points  of  symmetry.  With 
the  quadratic  form  cycle,  the  second  case  can  be  ignored  because  of  the  alternating  signs. 
However,  it  is  quite  possible  for  the  term  at  one  symmetry  point  to  be  merely  the  negative 
of  the  term  at  the  other  symmetry  point.  This  would  correspond  to  the  continued  fraction 
having  an  odd  period  and  there  would  be  a  symmetry  point  of  the  second  type  in  the 
continued  fraction  at  half-way.  However,  this  type  of  symmetry  does  not  generally  provide 
a  factorization  for  N. 

Lastly,  it  is  important  how  these  ambiguous  forms  fit  into  the  rest  of  the  class  group.  First, 
addressing  the  class  group  structure  requires  inverses.  Lemma  A. 2. 7  is  fairly  elementary  and 
is  probably  stated  somewhere  else.  Let  1  represent  the  form  in  the  principal  cycle  whose  first 
coefficient  is  1.  Let  F indicate  the  symmetric  reverse  of  F,  (a,b,  c)^1  =  (c,  b,  a).  Lemma 
A. 2. 7  justifies  this  notation: 

Lemma  A. 2. 7.  F  *  F~l  ~  1 
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Proof: 

Let  F  =  ax2  +  bxy  +  cy2.  Then  F ~1  =  cx 2  +  bxy  +  ay2.  Let  G  be  the  next  form  adjacent 
to  F~1,  that  is  G  =  ax2  +  b'xy  +  c'y2,  with  a  \  (b  +  b')  from  the  correspondence  with  continued 
fractions.  Composing  F  *  G,  n i  =  712  =  1  and  m  =  a,  so  that  the  first  coefficient  of  F  *  G  is 

1.  Therefore,  F  *  G  ~  1,  but  F ~  G,  so  F  *  F1-1  ~  1.  QED. 

Note  that  this  implies  that  the  square  of  a  symmetry  point  is  1. 

Theorem  A. 2. 8  was  probably  known  by  Shanks,  since  SQLIFOF  depends  highly  on  it,  but 
it  does  not  seem  that  he  states  this  explicitly  anywhere. 

Theorem  A. 2. 8.  An  equivalence  class  has  order  2  or  1  in  the  class  group  if  and  only  if  it 
is  ambiguous. 

Proof: 

Let  A  be  an  ambiguous  class.  Let  F  E  A.  Then  F  ~  F1-1,  so  that  F  *  F  ~  F  *  F1”1  ~  1. 
Therefore  F  *  F  is  in  the  principal  cycle,  so  that  A  has  order  2  or  1  in  the  class  group. 

Conversely,  assume  that  an  equivalence  class  A  has  order  2  or  1  in  the  class  group.  Let 

F  E  A.  Then  F*F  is  in  the  principal  ideal,  so  that  F*F  ~  (F*F)~1.  But  from  composition, 
it  is  clear  that  (F  *  F1)-1  ~  F~l  *  F~l.  So  F  *  F  ~  F~l  *  F~l.  Since  the  class  group  is 
associative,  composing  on  the  right  with  F  maintains  equivalence.  Therefore: 

(F1  *  F )  *  F  ~  (F  ^  *  F1-1)  *  F 
1  *  F  ~  F~1  *  (F1-1  *  F) 

F  ~  F1'1 

Therefore,  A  is  ambiguous.  QED. 

Certainly  the  class  group  structure  is  interesting,  but  it  is  now  possible  to  return  to  the 
problem  from  the  Morrison-Brillhart  algorithm  of  Example  6  with  =  3,  so  Qq  =  9  doesn’t 
provide  a  nontrivial  factor  of  N.  The  quick  explanation  is  that  if  you  square  the  quadratic 
form  with  first  coefficient  Q 3,  you  obtain  the  quadratic  form  with  first  coefficient  Qq.  Since 
the  principal  cycle  is  closed  under  composition,  it  seems  as  though,  and  perhaps  would  be 
convenient  if,  the  forms  in  the  principal  cycle  formed  a  group.  However,  the  problem  of 
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reduction  prevents  this: 

Example  14.  Consider  the  quadratic  form  F  =  (36,  70,  —3),  with  determinant  4- 1333.  Com¬ 
pare  (F  *  F)  *  F,  with  F  *  F  *  F,  where  the  difference  is  that  in  the  first  the  result  is 
reduced  after  the  first  composition.  F  *  F  =  (324,  —38,  —3)  and  the  very  next  adjacent  form 
(—3,68,59),  is  reduced.  F*  (—3,68,59)  =  (—12,70,9),  which  is  already  reduced.  However, 
without  reduction  F*F*F  =  F*  (324,  —38,  —3)  =  (729, 448,  —348).  When  this  is  reduced, 
the  first  reduced  form  found,  after  2  steps,  is  (9,56,  —61). 

Therefore,  the  principal  cycle,  with  the  operation  being  composition  followed  by  reduc¬ 
tion,  doesn’t  even  meet  the  requirements  for  being  power  associative.  However,  the  obser¬ 
vation  that  the  two  results  are  adjacent  forms,  and  that  the  second  reduction  took  one  step 
longer,  prompts  us  to  dig  a  little  deeper. 

Understanding  this  requires  what  Shanks  referred  to  as  infrastructure  distance.  For 
m  <  n,  and  for  xt  G  T,  the  terms  in  the  continued  fraction  in  (A. 5),  Shanks  defined 
infrastructure  distance  by 


Dw(xm,xn)  =  log(  xk)  (A. 13) 

k=m-\- 1 

Lenstra  [12]  adds  a  term  of  |  log (Qn/Qm)  to  this,  with  the  effect  that  the  resulting 
formulas  are  slightly  simplified  but  the  proofs  are  more  complicated  and  less  intuitive.  This 
definition  is  used  by  Williams  in  [26]. 

Since  the  quadratic  forms  are  cyclic,  in  order  for  the  distance  between  two  forms  to  be 
measured  consistently,  it  must  be  considered  modulo  the  distance  around  the  principal  cycle. 

Definition  A. 2.9.  Let  n  be  the  period  of  the  principal  cycle.  The  regulator  R  of  the  class 
group  is  the  distance  around  the  principal  cycle,  that  is, 

R  =  Df(F0,  Fn)  =  D{  1,  Fw) 


Therefore,  distance  must  be  considered  modulo  R,  so  that  Dp  is  a  map  from  pairs  of 
forms  to  the  interval  [0,  R)  G  M.  The  addition  of  two  distances  must  be  reduced  modulo  R 
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as  necessary. 

Further  analysis  of  this  distance  will  require  two  more  tools:  ideals  and  lattices.  In  order 
to  relate  to  continued  fractions,  the  ideals  will  be  in  Z [y/N]  =  {a  +  by/N  :  a,  b  G  Z}  and  the 
lattices  will  be  in  Q (y/N)  =  {a  +  by/N  :  a,  b  G  Q}  where  iV  is  a  non-square  positive  integer. 

A.  3  Ideals 

Remark:  The  ideals  in  Z [v/iV]  typically  correspond  only  to  quadratic  forms  of  discriminant 
4 N.  Note  that  if  N  =  1  (mod  4),  then  7L(y/N]  is  not  the  ring  of  integers  for  Q (y/N).  For 
N  =  1  (mod  4),  an  analysis  of  ideals  in  Z[v/^+1]  is  also  interesting,  but  will  be  avoided  in 
the  interest  of  simplicity.  Quadratic  forms  of  discriminant  N  =  1  (mod  4)  may  be  related 
to  ideals  in  Z [y/N]  via  first  multiplying  by  2  to  obtain  quadratic  forms  of  discriminant  AN. 
For  £  G  Q (y/N),  let  £  refer  to  the  conjugate  of  £  (i.e.  1  +  y/3  =  1  —  y/3). 

The  norm  of  a  number  in  Q (y/N)  is  A/"(£)  =  ££  G  Q. 

To  simplify  notation,  the  symbols  H,  /,  J,  and  K  will  consistently  be  ideals,  u  and  v 
will  be  elements  of  ideals,  a  and  f3  will  be  elements  of  Z [y/N],  £  and  £  will  be  elements  of 
Q (y/N),  and  C  will  be  a  lattice. 

Our  definition  of  an  ideal  is  the  same  as  in  any  other  commutative  ring  with  identity: 

Definition  A. 3.1.  A  subset  J  of  a  ring  R  is  an  ideal  if  for  u,v  e  I,  u±v  e  /  and  for  a  G  R, 
u-  a  G  /  (that  is,  /  is  closed  under  addition  and  multiplication  by  an  element  of  R).  If  /  is 
an  ideal  of  R  —  Z [y/N],  define  L(I)  to  be  the  least  positive  rational  integer  in  I. 

Describing  ideals  will  require  the  notation  for  the  lattice  generated  by  a  set.  If 

ai,a2,...,ak  G  Z [Vn], 

denote7  the  lattice  generated  by  these  as 

'Observe  the  difference  between  the  use  of  [...]  here  and  in  §A.l.  This  expression  is  completely  unrelated 
to  rational  approximations. 
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k 

[«1,  a2, oik\  =  {^2  niai  '■  ni  e  z}  (A.  14) 

i= 1 

Lemma  A. 3. 2  identifies  necessary  and  sufficient  conditions  for  a  set  in  h[y/N]  to  be  an  ideal. 

Lemma  A. 3. 2.  For  Q,  s,  A,  P  e  Z,  A  non-square  and  positive,  [Q,  sVN  +  P]  is  an  ideal  of 
the  ring  Z[\/A]  if  and  only  if  sQ  \  A f{sy/N  +  P ),  s  \  Q,  and  s  \  P. 

Proof:  Assume  that  /  =  [Q,s\fN  +  P]  is  an  ideal  of  Z [y/N].  Then,  choosing  a  = 
P  —  sVN,  Af(sy/N  +  P)  G  I.  Since  this  is  an  integer,  Q  \  A f(sy/N  +  P).  Choosing 
ol  =  \[N,  QVN  e  I .  Therefore,  s  \  Q.  Since  Q  \  A f(s\/N  +  P),  this  also  implies  that  s  \  P. 
Therefore,  a  could  have  been  chosen  a  =  P/s  —  \/N  so  that  Q  \  A f{s\fN  +  P)/s ,  so  that 
sQ  |  A f(sy/N  +  P). 

Conversely,  let  I  =  [Q,  sa/A  +  P]  and  assume  that  sQ  \  Af(s\fN  +  P),  s  \  Q,  and  s  \  P. 
Closure  under  addition  is  trivial.  To  see  that  /  is  closed  under  multiplication  by  an  element 
of  Z  [Vn\,  one  need  only  consider  multiplication  by  1  and  \/N,  since  they  form  a  basis  for 
Z[ VA]-  Multiplication  by  1  is  trivial.  For  \fN , 

qVn  =  —{sVn  +  P)  -  -Q 
s  s 

and  Q/s  and  —P/s  are  integers.  Also, 

P  —P2  +  s2N 

(sVN  +  P)Vn  =  sN  +  pVn  =  -(sVn  +  P)  +  ( - ^ - )Q 

s  sQ 

and  ^  and  (— are  integers.  QED 

If  s  =  1,  an  ideal  is  primitive.  Since  s  \  P  and  .s  |  Q,  ideals  that  are  not  primitive  will 
often  be  written  (s)[Q,  \fN  +  P],  Let  I  be  the  set  of  all  primitive  ideals. 

Represented  in  the  form  /  =  [Q,  yfN  +  P],  it  is  clear  that  \Q\  is  the  smallest  positive 
rational  integer  in  /.  Define 


L(I)  =  min{/nZ+} 


(A.15) 
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At  this  point,  it  is  possible  to  define  a  correspondence  between  quadratic  forms  (of  dis¬ 
criminant  A  =  0  (mod  4))  and  ideals  by: 

$vfl(F(x,  y)  =  Ax2  +  Bxy  +  Cy 2)  =  [A,  (^)2  -  AC  +  (A.16) 

p 2  _  AT 

<M[Q,  Vn  +  P])  =  F(x,y)  =  Qx 2  +  2 Pxy  +  (— —)y2  (A.16') 

and  define  a  reduced  ideal  as  an  ideal  corresponding  to  a  reduced  quadratic  form.  Note  that 
A  =  AN. 

For  example,  the  quadratic  form  (15,  2  •  12,  —1)  corresponds  to  the  ideal  [15,  a/159  +  12], 
The  one  potential  problem  that  immediately  becomes  apparent  is  that  while  [15,  a/159  +  12] 
and  [—15,  a/159  +  12]  are  the  same  ideal,  (15,2  •  12,-1)  and  (—15,2  ■  12,1)  are  different 
quadratic  forms.  However,  it  is  apparent  that  the  negative  sign  is  merely  carried  through 
composition  without  affecting  the  computations.  Since  each  of  these  forms  is  in  the  same 
location  within  its  respective  cycle,  this  difference  will  not  be  important  to  this  investigation 
of  composition  and  distance. 

and  may  be  combined  to  obtain 

-)  =  [Q,Vn-  p] 

and  $n,T  is  defined  in  the  related  obvious  way. 

If  A  =  [cq]  and  B  =  [$],  i  =  1,  2,  3...,  d,  then  it  is  clear  that  A  —  B  if  and  only  if  there 
exists  a  d  x  d  matrix  M  with  determinant  ±1  such  that: 


(a/  =  M(/3i) 


where  (a/  and  (/3j)  are  vectors. 

For  these  purposes,  the  most  important  operation  with  ideals  is  their  multiplication. 


Multiplication  is  defined  by 
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[<+]  '  [Pj\  =  iaiPj) 

For  example, 

/  =  [15,  \/l59  +  12]  •  [10,  Vl59  +  13]  =  [150,  lOv7!^  +  120, 15^159  +  195,  315  +  25^159] 

The  4th  component  is  the  sum  of  the  2nd  and  3rd,  so  it  is  unnecessary  for  describing  the 
ideal. 

Applying  the  matrix 

1  0  0 
0  1  0 
0  -1  1 

which  has  determinant  1,  subtracts  the  2nd  component  from  the  third  to  obtain 

I  =  [150, 10\/l59  +  120,  5a/159  +  75] 

The  matrix,  with  determinant  1, 

1  0  0 
0  1-2 
0  0  1 

will  subtract  twice  the  3rd  component  from  the  2nd  to  obtain 

I  =  [150,-30,5^159  +  75] 

Here  the  1st  component  is  a  multiple  of  the  2nd  and  is  thus  unnecessary.  The  answer  is 


simplified  to  obtain 
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I  =  [30,  5\/l59  +  75]  =  5[6,  ^159  +  15] 

The  process  of  multiplying  ideals  can  be  greatly  simplified  by  several  well-known  formu¬ 
lae8  [14]. 

Theorem  A. 3. 3.  Let  I  =  [Q,  \/N  +  P]  and  J  =  [Q',  \J TV  +  P'}  be  ideals  of  Qy/N.  Let 
C  =  a  =  If  gcd(Q,  P,  C)  =  gcd  {QfPfC)  =  1,  then! -J  =  s[q,VN+p], 

where 


s  =  gcd(Q,Q',P  +  P’) 

(A.17) 

h  =  gcd(Q,Q',C,  C',2) 

(A. 18) 

q  =  hQQ'/s 2 

(A. 19) 

p  =  P  (mod  Q/s) 

(A. 20) 

p  =  P'  (mod  Q'/ s) 

(A. 21) 

N  +  PP '  (mod  QQ'/s) 

(A. 22) 

Proof9: 

Consider  the  product: 

I  ■  J  =  IQQ',  QyfN  +  QPf  Q'y/N  +  Q'P ,  N  +  PP'  +  (P  +  P')Vn]  (A.23) 

The  smallest  integer  in  I  ■  J  may  be  found  by  considering  the  smallest  integers  that  may 
be  produced  taking  these  elements  pair-wise. 


L(I  ■  J )  =  gcd (QQ',  lcm(Q,  Q')(P 


lcm {Q,P  +  P’)Q’C  lcm (Q',P  +  P')QC 
b  p  +  p,  '  P  +  P' 


8In  [14],  (A. 22)  is  stated  as  (P  —  p)(P'  —  p)  =  n  +  tp  +  p2  (mod  QQ'/s),  but  in  this  case  t  =  0  and 
n  =  —N. 

9Some  of  the  arguments  were  taken  from  Buell’s  proof  in  [2]  concerning  composition  of  quadratic  forms. 
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Let  s  =  gcd(Q,  Q P  +  P'),  h  =  gcd(Q,  Q',  C,  C',  2),  w  —  hQQ' / s.  Let  f  ^  2  be  a  prime. 
Let  a,b,c,d,e,k  be  the  largest  possible  integers  such  that  /“  |  Q,  fb  \  Q\  fc  \  (P  +  P'), 
fd  \  C,  fe  \  C’,  fk  |  (P-P').  Then  fa+b  ||  QQ'-,  /maxM)+fc  ||  1  cm(Q,Q')(P  -  P'), 

fmax(a,c)+b+e—c  II  lcm (Q,P +P')Q'C'  1  fma x(b,c)+a+d—c  II  lcm(Q,,P +P')QC \ 

J  II  P+P'  j  ana  j  ||  P+P'  )' 

The  following  analysis  proves  that  if  /  ^  2,  the  maximum  exponent  of  /  in  L(I  ■  J ) 
is  a  +  6  —  min(a,  6,  c)  while  if  h  =  2,  then  the  maximum  exponent  of  2  in  L(J  ■  J)  is 
a  +  b  +  1  —  min(a,  6,  c),  while  if  h  =  1,  then  the  maximum  exponent  of  /  in  L(I  ■  J)  is 
a  +  b  —  min(a,6,  c).  As  this  is  broken  in  several  different  cases,  an  outline  of  the  proof  is 
helpful: 

1) a  =  0or6  =  0orc  =  0 

2)  a  /  0,6  7^  0,  and  c  /  0 
2-1)  f^2 

2.1.1)  a  +  d^b  +  e 

f\(p-n 

fup-n 

2.1.2)  a  +  d  =  b  +  e 

f\(p-n 

mp-n 

2.2)  /  =  2 

2.2.1)  a  +  d  b  +  e 

c  >  1,  k  >  1 
c  >  1,  k  <  1 
c  =  1 

2.2.2)  a  +  d  =  6  +  e 

c  >  1,  k  >  1 
c  =  1 
k  =  1 


Case  1)  If  a  =  0,  then  max(a,  c)+6  +  e  —  c  =  6  +  e>6,  max(6,  c)+a  +  d  —  c>6  and 
a  +  6  =  6,  so  the  maximum  exponent  for  /  in  L(I  ■  J)  is  6  =  a  +  6  —  min(a,  6,  c).  Similarly,  if 
6  =  0,  then  the  maximum  exponent  is  a  =  a  +  6  —  min(a ,  6,  c). 

Assume  c  =  0.  fa+d  ||  QC  =  N  —  P2  and  fb+e  ||  Q'C'  =  N  —  (P')2,  subtracting, 
/min (a+d,b+e)  j  (p2  _  (p/)2)  =  (p  +  p/)(p  _  p').  Since  c  =  0,  then  /min(“+<i’6+e)  |  (P  -  P'). 

Therefore,  Jmax(a>&)+mm(a+d>b+e)  |  lcm((5,  Q')(P  —  P').  However,  max(a,  6) +  min(a  +  d,  6  +  e)  > 
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max(a,  b )  +  min(a,  b)  =  a  +  b.  Therefore,  the  maximum  exponent  for  /  in  L(I  ■  J)  is 

min(a  +  b ,  max(a,  c)  +  b  +  e  —  c,  max(6,  c)  +  a  +  d  —  c)  =  min(a  +  b,  a  +  b  +  e,  a  +  b  +  d) 

=a+b=a+b—  min(a,  b ,  c) 

Case  2.1.1)  Assume  c  ^  0,  /  ^  2,  and  a  +  d  ^  b  +  e.  Then,  /mm(a+d,fe+e)  ||  (P2  —  (P')2)  = 
(P  +  P')(P  -  P')-  If  /  |  (P  -  P'),  then  /  |  2P  and  /  |  2P'.  Since  /  ±  2,  this  gives  /  j  P, 
/  |  P'.  Then,  d  =  e  =  0.  Also,  c  <  min(a,  b).  Then,  the  maximum  exponent  for  /  in  L(I  ■  J) 

is 

min(a  +  b,  max(a,  b)  +  min(a,  b)  —  c,  max(a,  c)  +  b  —  c,  ma x(b,  c)  +  a  —  c)  = 

a  +  6  —  c  =  a  +  6  —  min(a,  b,  c ) 

If  /  j  (P  —  P')  then  c  =  min(a  +  d,  b  +  d)  and  the  maximum  exponent  for  /  in  L(J  •  J)  is 

min(a  +  b ,  max(a,  6),  max(a,  c)  +  b  +  e  —  c,  max(b,  c)  +  a  +  d  —  c) 

If  a  =  min(a,  6,  c),  then  this  is  min(&,  c+b+e  —  c ,  max(6,  c)+a+d— c)  —  b  —  a+b— min(a,  b,  c ). 
The  case  is  similar  if  b  =  min(a,  6,  c).  If  c  =  min(a,  6,  c),  then  since  c  =  min(a  +  d,  b  +  e), 
this  gives  c  =  ruin  (a,  b).  Then  the  maximum  exponent  is 

min(max(a,  b),a  +  b  +  e  —  c,  b  +  a  +  d  —  c)  =  max(a,  b)  =  a  +  b  —  min(a,  b)  =  a  +  b  —  min(a,  6,  c) 

Case  2.1.2)  Assume  c  ^  0,  /  /  2,  but  a  +  d  —  b  +  e.  As  before,  if  /  ]  (P  —  P'), 
then  d  =  e  =  0.  In  this  case  also  a  =  b.  Assume  c  <  a.  Note  that  fa~c  \  (P  —  P')  and 
max(a,  b)  +  min(a  +  d,  b  +  e)  —  c  =  a  +  b  —  c.  Then  the  maximum  exponent  is 


min(a  +  b ,  max(a,  c)  +  b  —  c,  max(6,  c)  +  a  —  c)  —  a  +  b  —  c  =  a  +  b  —  min(a,  6,  c). 
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Alternately,  assume  c  >  a.  Then  for  some  k,  fk  ||  (P  —  P').  The  maximum  exponent  is 
min(a  +  6,  max(a,  b )  +  k,  6,  max (6,  c)  +  a  —  c)  =  b  =  a  +  b  —  a  =  a  +  b  —  min(a,  b ,  c ) 
Conversely,  assume  f  \  (P  —  P').  Then  c>a  +  d,  =  b  +  e  and  the  maximum  exponent  is 

min(a  +  b ,  max(a,  b),  max(a,  c)  +  b  +  e  —  c,  max(6,  c)  +  a  +  d  —  c)  =  min(max(a,  b),  a  +  d) 

=  max(a,  b)  =  a  +  b  —  min(a,  b)  =  a  +  b  —  min(a,  6,  c) 

Case  2.2.1)  Let  /  =  2.  Assume  a  +  d  ^  b  +  e.  Then  2mm(c’fc)  ||  2P  and  2mm(c,fc)  ||  2P',  so 
that  2min(c’fcC1  ||  Pi  p'm  If  c  >  1  and  k  >  1,  then  d  =  e  =  0  and  as  before  the  largest  exponent 
is  a  +  b  —  min(a,  b,  c ).  Assume  c  >  1,  k  <  1.  Then  k  =  1  and  c  =  min(a  +  d,  6  +  e)  —  1.  The 
largest  exponent  is  then 

min(a  +  b,  max(a,  b)  +  1,  max(a,  c)  +  b  +  e  —  c,  max(6,  c)  +  a  +  d  —  c) 

If  a  <  min(6,  c)  this  reduces  to  &  +  min(e,  1)  =  a  +  6  +  min(e,  1)  —  min(a,  b,  c).  c+l<a  +  d< 
c  +  d,  so  d  >  1.  Note  that  if  e  >  1  then  this  is  a  special  case  and  h  =  2.  If  e  =  0,  it  is  the 
same  as  before.  The  cases  when  b  <  min(a,  c)  are  similar. 

If  c  <  min(a,  b ),  this  exponent  is  min(max(a,  6)  +  l,  a  +  b+e  —  c,  b+a  +  d  —  c ).  Without  loss 
of  generality,  assume  a  +  d>&  +  eso  that  c  =  a  +  d  —  1  >  c  +  d—  1,  so  that  d  =  1  and  a  =  c. 
Then  the  exponent  is  min(6+l,  6+e,  b+d)  =  6+min(l,  e,  d)  =  a+6+min(l,  e,  d)— min(a,  b,  c ). 
Note  again  that  if  e  >  1  and  d>  1,  then  this  is  the  special  case  where  h  =  2.  Otherwise,  it 
is  the  same  as  before. 

Assume  c  =  1.  Then  k  >  1.  Then  the  exponent  is 


min(a  +  6,  max(a,  b)  +  min(a  +  d,  b  +  e)  —  1,  a  +  b  +  e  —  1,  a  +  b  +  d  —  1). 
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If  d  =  0  or  e  =  0,  h  =  1  and  this  is  a  +  6  —  1  =  a  +  6  —  min(a,  6,  c).  Otherwise,  h  =  2  and 
the  exponent  isa  +  6  =  a  +  6  +  l  —  min(a,  6,  c ). 

Case  2.2.2)  Lastly,  assume  that  c  /  0  but  a  +  d  —  b  +  e.  For  some  k,  2k  ||  (P  —  P'). 
c  +  k  >  a  +  d.  If  c  >  1  and  k  >  1,  then  d  =  e  =  0,  a  =  b.  The  exponent  is  then 
min(2a,  a  +  k,  max(a,  c)  +  b  —  c).  If  c  >  a,  this  is  min(2a,  a  +  k,  a)  =  a  =  a  +  b  —  min(a,  b,  c ). 
If  c  <  a,  the  exponent  is  min(2a,  a  +  k,2a  —  c)  =  2a  —  c  =  a  +  b  —  min(a,  b,  c ). 

Alternately,  if  c  =  1,  then  k  >  a  +  d  —  1  and  the  exponent  is 

min  (a  +  6,  max(a,  b)  +  k,a  +  b  +  e  —  l,b  +  a  +  d  —  1)  =  min(a  +  b,  a  +  b  +  e  —  l,a  +  b  +  d—  1). 

If  e  >  0  and  d  >  0,  h  =  2  and  this  exponent  isa  +  6  =  a  +  6  +  l  —  min(a,  b,  c).  If  e  =  0  or 
d  —  0,  h  —  1  and  this  isa  +  6—  l  =  a  +  6  —  min(a,  b,  c ). 

If  k  —  1  then  c  >  1  and  specifically  c  >  a  +  d  —  1  =  b  +  e  —  1.  If  c  <  min(a,  6),  then 
d  =  e  =  1  so  that  h  =  2,  c  =  a  =  b  and  the  exponent  is 

min(2a,  a  +  l)  =  a  +  l  =  a  +  5  +  l  —  min(a,  6,  c) 

If  a  <  min (b,  c),  then  d  >  e  and  the  exponent  is 

min(a  +  b,  b  +  1,  b  +  e,  max(6,  c)  +  a  +  d  —  c)  =  min(6  +  1,  b  +  e) 

If  e  >  1,  then  d  >  e  >  1,  so  h  =  2  and  in  this  case  the  exponent  is  6+1  =  a+6+1— min(a,  6,  c). 
If  e  =  0,  h  —  1  and  in  this  case  the  exponent  is  b  =  a  +  6  —  min(a,  6,  c). 

Therefore, 

L{I  ■  J)  —  hQQ'/s 

so  that  for  /  ^  2,  the  highest  exponent  is  a  +  b  —  min(a,  6,  c)  and  for  /  =  2,  the  highest 
exponent  is  a  +  b  —  min(a,  6,  c)  if  h  =  1  and  a  +  6  +  1  —  min(a,  6,  c)  if  h  =  2.  Note  that  this 
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is  still  divisible  by  s.  After  that  s  is  factored  out,  the  result  is  q  —  hQQ' / s2 . 

There  are  integers  t,  u,  and  v  such  that  tQ  +  uQ'  +  v(P  +  P')  =  s.  First,  consider 
divisibility  by  s.  This  is  trivial  for  every  term  except  N  +  PP' .  By  the  definition  of  s, 

P  +  P'  =  0  (mod  s) 

P'  =  —  P  (mod  s) 

PP '  =  -P2  (mod  s) 

N  +  PP1  =  N  —  P2  (mod  s) 

and  since  s  \  Q  and  Q  \  (N  —  P2),  s  \  N  +  PP'. 

The  linear  combination  of  the  last  three  elements  with  coefficients  t,  u,  and  v  respectively 
is: 


sVN  +  tQP '  +  uQ’P  +  v(N  +  PP') 

so  that  it  is  evident  that  after  s  is  factored  out,  the  remaining  ideal  is  primitive.  Since  this  is 
the  element  of  I  ■  J  with  the  smallest  coefficient  of  y/N,  clearly  p  =  t(Q/s)P'  +  u(Q'/s)P  + 
v(N  +  PP')/s,  modulo  L(I  ■  J).  Then, 

p  =  t(Q/s)P'  +  (s  —  tQ  —  v(P  +  P'))P/s  +  v(N  +  PP')/s 
=  P  +  v(N  —  P2)/s  (mod  Q/s) 

=  P  (mod  Q/s ) 

since  Q  \  (N  —  P2).  By  symmetric  arguments,  p  =  P'  (mod  Q' / s). 


To  prove  (A. 22),  consider: 
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(P  +  P')sp  =  (P  +  P')(tQP'  +  uQ'P  +  v(N  +  PP') 

=  (P  +  P')(tQP '  +  uQ'P)  +  (P  +  P'){N  +  PP> 

=  (P  +  P')(tQP'  +  uQ'P)  +  (s  -  tQ  -  uQ')(N  +  PP') 

=  s(JV  +  PP')  +  fQ((P')2  -  P)  +  mQ'(P2  -  N) 

=  s(N  +  PP')  (mod  QQ') 

Therefore,  (P  +  P')p  =  N  +  PP'  (mod  QQ'/s).  QED 
Observe  that  when  h  =  1  (A.  19)  could  be  restated  as 

L(I  ■  J)  =  L{I)L(J)/s2  (A. 24) 

remembering  that  L(I)  is  defined  as  the  smallest  positive  rational  integer  in  I.  This  equation 
is  proven  in  [26]  and  will  be  useful  later. 

Also  observe  that  for  h  =  1,  the  equations  describing  the  product  of  two  ideals  correspond 
exactly  to  the  composition  of  two  quadratic  forms.  Shanks  notes  this  in  [21].  Therefore,  the 
equations  concerning  distance  and  multiplication  of  ideals  will  correspond  to  distance  and 
composition  of  quadratic  forms. 

The  case  when  h  —  2  connects  composition  of  quadratic  forms  of  discriminant  =  1 
(mod  4)  to  multiplication  of  ideals.  If  P  and  G  are  two  quadratic  forms  with  discriminant 
N  =  1  (mod  4),  then  2 P  and  2 G  have  discriminant  4A^  and  correspond  to  ideals  I^f  and 
J2g  in  h[y/N].  Multiplying,  h  —  2  and  Pf  •  he  =  P(f*g)-  Therefore,  although  this  case  will 
not  be  considered  further,  it  is  readily  seen  that  the  distance  formulas  derived  from  ideals  in 
Z[y/ff\  will  still  correspond  to  composition  of  quadratic  forms  of  discriminant  =  1  (mod  4). 
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A. 4  Lattices 

Consider  lattices  in  Q (y/N)  x  Q (y/N).  Define  L  as  the  set  of  all  lattices  in  Q (y/N)  x  Q (\/N). 
Define  the  map  M  :  Q (y/N)  — ■»  Q (y/N)  x  Q (y/N)  by 

m(o  =  <£,£> 

Multiplication  in  Q (y/N)  x  Q (y/N)  is  defined  component-wise,  that  is,  (£,£')  ■  ((,(')  — 
(£C,£'C'>,  so  that  it  is  clear  that  M  is  homomorphic  and  one-to-one. 

Distance  will  relate  to  a  concept  called  a  minimum: 

Definition  A. 4.1.  For  a  vector  v  =  (vi,v2,  ...Vd),  the  normed  body  of  v,  TZ(v)  is  the  set 

TZ{v)  =  {{x-l,x2 ,...xd)  ■  Xi  G  M,|xj|  <  \vi\,i  =  1,2 
Abusing  notation,  denote  7£(£)  =  7Z((£,£)). 

A  number  ^  (or  actually  the  corresponding  vector)  is  a  minimum  of  C  if  H(£)  D  C  —  {0}, 
where  0  is  the  vector  (0,  0). 

A  lattice  C  is  reduced  if  1  e  C  and  1  is  a  minimum. 

For  this  case  with  d  —  2,  the  normed  body  is  a  rectangle  in  M2.  Note  that  for  £  G  Q (\/N) 
the  normed  body  H(£)  has  area  equal  to  four  times  the  absolute  value  of  the  norm  |A/"(£)|. 

To  avoid  unnecessary  generality,  this  investigation  will  focus  specifically  on  the  lattices 
corresponding  to  ideals.  Specifically,  for  the  primitive  ideal  /  =  [Q,\/N  +  P],  define  the 
associated  lattice  containing  1  in  Q(a /N)  as  £,  =  [i  ,(Vn  +  p)/q\. 

Conversely,  to  each  lattice  containing  1  in  Q(a /N)  there  is  an  associated  primitive  lattice 
(which  may  or  may  not  be  an  ideal)  in  Z [v/iV] .  Equation  (A. 15)  defined  the  function  L.  In 
a  similar  fashion,  for  a  lattice  £,  define 
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L(£)  =  min{n  e  Z+  :  n£  C  Z [ViV]}  (A. 25) 

Then  if  L(£)£  is  an  ideal  of  Z [y/N]  it  is  the  primitive  ideal  associated  to  a  lattice  £.  Note 
that  if  an  ideal  I  is  associated  to  a  lattice  £/,  then  L(I)  =  L(£j).  Define 

Vn  +  P})  =  [1,  (Vn  +  P)/Q } 

and 

$l,i(£)  =  L(£)£ 

Note  that  for  some  lattices  C,  may  not  actually  be  an  ideal.  Lemma  A. 4. 2  provides 

conditions  for  it  to  be  an  ideal  sufficient  for  this  analysis: 

Lemma  A. 4. 2.  Let  I  be  a  primitive  ideal  and  let  £  =  If  £'  is  a  lattice  with  basis 

{1,£}  and  for  some  9,  9£!  =  £,  then  J  =  is  a  primitive  ideal  and 

{L{I)9)J  =  (L(J))I 

Proof:  Let  /  =  [Q,  \/N  +  P],  Then  £  =  [1,  (\/N  +  P)/Q).  The  statement  that  9£!  =  £ 
requires  that 

1  1 
9  =  T  _ 

[(VN  +  P)/Q 

where  T  is  a  2  x  2  matrix  with  determinant  ±1.  Multiplying  by  L(I)  =  L(£)  =  Q  and 
L(J)  =  L(£')\ 

L(£')  Q 

Q9  =  L(£')T 

L{£')Z  ( VN  +  P ) 
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so  that  (. L(I)6)J  =  (L(J))/.  Therefore,  J  is  an  ideal.  It  is  primitive  by  the  definition  of 
^L,!-  QED 

For  an  example  of  minima,  consider  the  lattice  [1,  \/l59  — 12],  7^.(1)  is  a  square  with  sides 
of  length  2  centered  at  the  origin  and  a  simple  graph  demonstrates  that  0  is  the  only  point 
in  the  lattice  and  contained  in  this  square.  Therefore  1  is  a  minimum,  a/159  —  12  is  also  a 
minimum.  7\!.(a/159  +  12)  is  a  narrower  and  taller  rectangle  also  centered  at  the  origin. 

Given  two  minima,  it  is  important  to  be  able  to  determine  whether  or  not  there  is  another 
minimum  between  them.  In  vector  format,  if  (xi,yi)  and  (x2,  y2)  are  minima  with  /i|  >  |x2| 
and  \yi\  <  \ y2\,  these  two  minima  are  adjacent  if  there  does  not  exist  another  minima  (x^,  yf) 
such  that  /2|  <  | £3 1  <  |xi|  and  \yx\  <  \y3\  <  \y2\. 

Voronoi  developed  a  method  (and  a  theorem)  concerning  adjacent  minima  ([4],  [26]). 

Theorem  A. 4. 3.  Let  C  be  a  lattice  with  {£,£}  as  a  basis,  where  £,  £  €  Q (VN)  arid  suppose 
that  £  >  £  >  0.  Then  £  and  £  are  adjacent  minima  of  C  if  and  only  if  |£|  >  |£|  and  ££  <  0. 

Proof:  Assume  £  and  £  are  adjacent  minima.  Since  they  are  both  minima,  |£|  >  |£|,  or 
else  £  would  not  be  a  minima.  Also  0  <  £  —  £  <  (.  Since  (  is  a  minima,  this  requires  that 
|C  —  >  |C|-  If  C  and  C  had  the  same  sign,  this  would  not  be  possible.  Therefore,  <  0. 

Conversely,  assume  that  |£|  >  |C|  and  <  0.  Assume  that  £  is  not  a  minimum  of  C. 
Then  there  exists  some  u  G  Q (\/N)  such  that  /j  <  £  and  |oJ|  <  |£|.  Since  lu  =  a£  +  b(  for 
some  a,  6  6  Z,  |a£  +  &£|  <  £  and  |a£  +  5£|  <  |£|.  If  ab  =  0,  then  either  a  —  0  or  b  —  0.  If  a  =  0, 
then  the  second  statement  contradicts  the  hypothesis.  If  b  =  0,  then  the  first  statement  gives 
£  <  £,  clearly  false.  However,  if  ab  >  0  then  |a£  +  b(j\  >  £  and  if  ab  <  0,  then  since  ££  <  0, 
|a£  +  bf  |  >  |£|.  Therefore,  £  must  be  a  minima.  By  similar  reasoning,  £  must  be  a  minima. 

Concerning  adjacency,  assume  that  there  is  another  minima  cu  between  £  and  £.  Since 
00  =  a£  +  b(  for  some  a,  b  e  Z,  £  <  |a£  +  &£!<£  and  ICI  <  |a£  +  b(\  <  |£|.  Since  £  >  £  >  0, 
the  first  statement  requires  that  6  =  0  and  then  the  second  statement  simplifies  to  |a|  <  1, 
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requiring  that  a  =  0  and  providing  a  contradiction.  Therefore,  £  and  £  are  adjacent  minima. 

QED 

From  the  previous  example,  it  is  now  possible  to  check  that  £  =  1  and  £  =  y/159  +  12 
are  indeed  adjacent  minima. 

The  idea  that  will  actually  connect  to  continued  fractions  (and  distance)  is  the  search  for 
a  sequence  of  adjacent  minima.  This  sequence  is  formed  by  relating  different  lattices.  The 
following  Lemmas  are  due  to  Williams  [26]. 

Lemma  A. 4. 4.  Let  C  and  £  be  reduced  lattices.  If  ££'  =  C.  then  £  is  a  minimum  of  C. 

Proof:  Since  1  G  £ ,  £  G  C.  If  £  is  not  a  minimum  of  £,  then  there  exists  a  £  G  C  such 
that  £  ^  0  and  |£|  <  |£|  and  |£|  <  |£|.  Let  [3  =  £/£,  so  that  /3  G  £ .  \(3\  =  [£/£[  <  1  and 
\(3\  =  [£/£[  <  1,  contradicting  the  fact  that  C  is  reduced.  Therefore,  £  is  a  minimum  of  C. 

QED 

Now  consider  the  converse  of  this  statement.  Note  that  |_^J  denotes  the  floor  of  x. 

Lemma  A. 4. 5.  Let  C  =  [1,£],  where  1  and  £  are  adjacent  minima  of  C  with  1  >  £  >  0.  Let 
£  =  (1  /£)£.  Then  £  is  a  reduced  lattice. 

Proof:  £  =  (1  /£)  [1,  £]  =  [1  /£,  1]  =  [l/£  —  |_1/£J ,  1] ,  so  that  1  G  £.  It  is  sufficient  to 
show  that  1  and  £'  =  l/£  —  |_1/£J  are  adjacent  minima.  First,  1  and  £'  are  a  basis  for  £ 
and  1  >  £'  >  0.  Since  0  <  £  <  1,  |_1/£J  >  1.  Since  £  <  0,  £'  =  l/£  —  [_1  / £J  <  0  —  1  =  — 1. 
Thereby  satisfying  both  the  requirement  that  £  •  1  <  0  and  the  requirement  that  |£  |  >  1. 
Therefore,  by  Theorem  A. 4. 3,  1  and  £'  are  adjacent  minima  of  £  and  thus  £  is  a  reduced 
lattice.  QED 

Actually,  these  proofs  provide  a  bit  more  by  actually  finding  the  minimum  adjacent  to  1 
in  the  new  lattice.  The  next  Lemma  makes  use  of  this  minimum  [26]: 

Lemma  A. 4. 6.  Let  C,  £ ,  £,  and  £'  be  as  above.  Let  £  be  the  minimum  adjacent  to  £  other 
than  1  in  C.  Then  £  =  ££'. 


87 


Proof:  ££'  =  £(l/£-|_l/£J)  =  !— ^L!/CJ ,  so  that  [£,££']  is  a  basis  for  C.  Since  1  >  £'  >  0, 
£>££'>  0.  Since  |£'|  >  1,  |££'|  >  |£j.  Since  £'  <  0,  £  •  ££'  =  (£)2£'  <  0.  Therefore,  by 
Theorem  A. 4. 3,  £  and  ££'  are  adjacent  minima.  Since  ££'  ^  1,  £  =  ££'.  QED 

Observe  that  by  a  similar  process,  one  could  find  a  reduced  lattice  C"  =  !/£'£',  etc.  Then 
C"  =  1  /(££')£.  To  generalize,  define  £  =  £1  and  C  —  C\  and  this  is  a  sequence  of  reduced 
lattices  and  their  minima.  A  chain  of  adjacent  minima  of  C\  may  be  defined  by 


n—  1 

0n  =  n& 


and  then 


(A. 26) 


erXn  =  Cl  (A. 27) 

Since  each  is  a  reduced  lattice,  by  Lemma  A. 4. 4  each  9n  is  a  minimum  of  C\. 
Although  it  is  not  true  in  higher  dimensions,  it  is  fairly  trivial  in  2-d  that  this  chain  of 
adjacent  minima  provides  a  complete  (although  infinite)  list  of  the  minima  with  ^-coordinate 
between  0  and  1. 

Lemma  A. 4. 7.  Let  (0,  0)  be  a  minimum  of  a  lattice  C,  with  0  <  0  <  1.  Then  for  some  n, 
0  =  0n,  where  0n  is  defined  by  equation  (A. 26) 

Define  distance  in  terms  of  this  chain  of  minima  by 


Dh(Cn,Cm)  =  log(0n/0m)  (A. 28) 

It  will  become  readily  apparent  that  the  subscript  L  is  unnecessary,  but  it  provides  clarity 
for  now.  Before  continuing  it  is  appropriate  to  provide  an  example  of  these  concepts.  First, 
as  a  reference,  consider  the  steps  for  the  continued  fraction  expansion  of  Vl59  —  12  and  the 
quadratic  form  distances  D f  covered  to  the  end  of  each  step: 


_  1  ^159+12  _  i  yT59-3 

1  VT59-12  15  15  ’ 

_  15  _  \/T59+3  _  i  ,  \/l59— 7 

2  x/159— 3  10  “r  10  ’ 

_  10  _  VI59+7  _  i  ,  \/l59— 4 

3  V159-7  11  11  ’ 

_  11  _  v/159+4  _  i  ,  \/l59— 9 

X4  “  \/l59— 4  “  13  _  1  "T"  13  ’ 

_  13  _  CI59+9  _  q  ,  V159-9 

X5  ~  x/159— 9  “  6  “  °  6  ’ 


JDF(F0,F1)  =  log(^M) 

JDF(F0,F2)  =  log(^fP) 
Dw(F0,F3)  =  log(^M±25) 
JDF(Fo,F4)  =  log(^f|±^) 
Dw(Fq,F5)  =  iog(5y2M+63)_ 


The  continued  fraction  corresponds  to  quadratic  forms  which  correspond  to  ideals,  which 


are  associated  with  lattices  that  contain  1.  In  this  case,  the  lattice  associated  with  x\  is 
C\  =  [1,  1/aq]  =  [1,  \/l59  —  12]  and  l/x\  is  a  minimum  adjacent  to  1  in  C\.  From  here: 


^2  “  vA-12^1  “  1^155-12’ -1] 

=  Tifc^2  =  i^lS-3’1] 


[nsti?,  i]  =  [i,  u|M]  =  [i,  i/x2] 


and  it  is  apparent  that  this  same  pattern  of  correspondence  will  continue,  that  is 


—  <£t,i(<£]i,]l(T7i))  —  [1,1  /xn]  —  Cn  (A. 29) 

from  which  it  is  also  apparent  that  the  sequences  of  lattices  will  be  periodic. 

Computing  equation  (A. 26),  for  example,  93  =  (\/l59  —  12)(  v'^|~-)  =  13  —  vT59.  With 
0i  =  l, 

D{CU  £3)  =  log(l/(13  -  Vl59))  =  log((Vl59  +  13)/10) 

It  is  readily  apparent  that  the  definition  of  distances  in  lattices  corresponds  to  the  defi¬ 


nition  given  for  quadratic  forms.  Note  that  these  distances  must  still  be  considered  modulo 
R ,  the  regulator,  since  the  sequence  of  lattices  is  still  cyclic. 
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A. 5  The  Generalized  Distance  Formula 

Going  back  to  ideals,  note  that  if  I\  =  L(Ci)Ci  and  In  =  L(Cn)Cn  is  another  ideal  corre¬ 
sponding  to  a  lattice  later  in  the  same  sequence,  then 


enc 


n 


Cl 


L(£1)L(£r),)^£„  =  L(Ci)L(Cn)Ci 

(L(A)0n)/n  =  {L(Cn)h 


(A. 30) 


where  once  again,  the  distance  (this  time  between  ideals)  is  given  by  D(Ii,In)  =  —  log(0n). 
Now,  this  definition  of  distance  is  well  and  good  for  reduced  ideals,  but  as  of  yet,  it  hasn’t 
been  applied  it  to  non-reduced  ideals.  To  relate  the  definitions  of  reduced  lattices  and 
continued  fractions  observe  that  the  definition  of  a  reduced  continued  fraction  implies  that 
for  a  term  Xi  =  ,  being  reduced  equates  to 


Vn  +  Pj-i 

Qi 


>  l 


Vn  -  Pi-i 

o  <  - —  ?  1  <  1 


so  that  it  is  clear  that  if  Cx  =  [1,1/x],  then  1  and  x  are  adjacent  minima  and  the  lattice 
is  reduced.  The  process  of  dealing  with  a  non-reduced  lattice  correlates  to  the  process  of 
reducing  a  continued  fraction  as  demonstrated  in  the  proof  of  Lemma  A. 1.8.  See  [26]  for  a 
more  general  Lemma. 


Lemma  A. 5.1.  Let  I  be  any  primitive  ideal  in  h[y/N].  There  exists  a  reduced  ideal  In  and 
a  9n  G  /  such  that 
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(■ L{I)Qn)In  =  ( L(In))I  (A. 31) 

Proof: 

Let  I  =  [Q,y/N  +  P\.  Then  the  associated  lattice  is  Cj  =  [1,  AA+P]  =  [l,£i].  If  / 
is  reduced,  In  —  /,  u  =  L(J),  and  the  proof  is  done.  If  /  is  not  reduced,  then  Cj  is  not 
reduced.  Without  loss  of  generality,  assume  that  0  <  <  1  (since  otherwise  it  would  just 

have  to  be  reduced  by  an  integer.).  Let  C2  =  1/£l£j  =  [l/£,  1]  =  [l,l/£  —  |_l/£  —  1/2J ] . 
Then  £,\£2  =  Cj.  Continuing  in  similar  manner10,  by  Lemma  A. 1.8  and  the  correspondence 
between  lattices  and  continued  fractions  for  some  n,£n  reduced,  and  thus  Cn  reduced.  As  in 
(A. 26),  set 


so  that 


n—  1 

0n  =  n& 
2—1 


6nCr 


£ i 


Then  ( L{I)Qn)In  =  ( L(In)I .  QED 

Let  A,  J\  be  reduced  primitive  ideals.  Let  K\  be  the  primitive  ideal  found  by  multiplying 
I\  and  -J\  and  removing  a  factor  and  let  s  be  the  factor  removed,  so  that  (s)Ki  =  I\Ji,  s  G  Z. 
By  Lemma  A. 5.1  there  exists  a  reduced  ideal  Kj  and  a  A j  G  K\  such  that 

{L{K1)\j)Kj  =  (L(A^))Ad  (A. 32) 

corresponding  to  D(Ki,Kj)  =  «*log(Aj). 

Let  In  ~  / 1  and  Jm  ~  J\  and  let  Hi  be  the  primitive  ideal  found  by  multiplying  In  and 

10Note  that  it  is  irrelevant  whether  or  not  the  second  components  of  the  intermediate  lattices  are  either 
minima  or  adjacent  to  1.  Also  note  that,  as  in  Lemma  A. 1.8,  the  formula  would  change  slightly  when  the 
denominators  get  small. 
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Jm  and  removing  a  factor  and  let  t  be  the  factor  removed,  so  that  (t)Hy  = 
Lemma  A. 5.1  there  exists  a  reduced  ideal  Hk  and  a  qk  E  K\  such  that 

(L(Hi)r]k)Hk  =  (L<yHk))Hi 

corresponding  to  D(Hi,  Hk)  =  —  \og('qk). 

Also,  there  exist  minima  /in  and  (pm  in  the  lattices  corresponding  to  I\  and  J\ 
such  that 


(■ L{h)nn)In  =  ( L(In))h 


and 


corresponding  to  D(Ii,In)  =  —  log(/in)  and  D(Jk,  Jm )  =  —\og(4>n 
By  combining  (A. 24)  and  (A.32)-(A.35): 


(W)ifj  =  (Wwr)  K> 

_  ( L{Hk)L(K.j)\  ,  j 
V  HK^XjS  J  llJl 


V  L(Ki)\jsL(In)L(Jm)  I  ±n 


In.Jrr 


_  {  L(H Sflntfin 

-  \  XjL(In)L(Jm)  I  ±nUm 

__  f  L(Hk)L(Kj)sfin<f>  mt 
~  ^  AjL(/„)L(Jm) 


‘  Xjtnih) 

_  f  L(Kj)skln<l>m.'nk 


In  'I : 

Hi 

Bx 


Hk 


s  j  I  n  (.)  nl  Ijk 
t\j 


m,  te  Z.  By 


(A. 33) 


,  respectively, 


(A. 34) 


(A. 35) 


Set 
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and  then 

(L{Kj)^)Hk  =  (L(Hk))Kj  (A.36) 

Since  Kj  and  Hk  are  reduced,  by  Lemma  A. 4. 4  if  is  a  minimum  of  the  lattice  C k:i  ,  so 
that  for  some  n,  if  =  9n.  Therefore, 

D(Kj,  Hk)  =  -  log(^)  =  -  log (nn)  -  log (<f>m)  -  log (rjk)  +  log(Aj)  -  log(s/t ) 

—  D(I\,  In )  +  D(  Ji,  Jm)  +  C 

where  C  =  D(Hi,  Hj)  —  D(Ki,  Kj)  +  logit/ s)  will  be  small  compared  to  D(Kj,  Hk)  for  m,n 
large. 

By  the  correspondence  between  multiplication  of  ideals  and  composition  of  quadratic 
forms,  this  result  may  be  restated  in  terms  of  forms: 

Theorem  A. 5. 2.  If  Fx  ~  Fn  are  equivalent  forms  and  G i  ~  Gm  are  equivalent  forms  and 
Dp  i  is  the  reduction  distance  for  F\  *  G\  and  Dp$  is  the  reduction  distance  for  Fn  *  Grn  and 
s  and  t  are  the  factors  cancelled  in  each  respective  composition,  then 


D(F i  *  GuFn  *  Gm)  =  D{FU  Fn)  +  D{GU  Gn)  +  C 
where  (  =  DPt2  -  Dpj  +  log(t/s). 

Example  6  from  the  Morrison-Brillhart  algorithm  is  in  the  principal  cycle.  By  Theorem 
A. 5. 2  when  a  form  F  is  composed  with  itself,  the  distance  from  1  to  F  is  roughly  doubled, 
d{l,F2)  =  2d(l,F)  +  (.  Therefore,  the  index  is  roughly  doubled,  since  distance  is  roughly 
proportional  to  the  difference  in  indices,  so  that  F3  *  F3  —  F6,  and  Qq  =  9  is  the  square  of 
Qs  =  3. 

Since  the  square  of  any  symmetry  point  has  first  coefficient  1,  observe  that  if  the  distance 
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around  some  cycle  were  unrelated  to  the  distance  around  the  principal  cycle,  then  this 
result  would  be  affected  by  which  symmetry  point  this  distance  was  referenced  from.  From 
Definition  A. 2. 9  R  =  D(F0,Fn )  in  the  principal  cycle.  At  this  point,  it  is  clear  that  the 
distance  in  other  cycles  must  be  the  same. 

Lemma  A. 5. 3.  Let  A  be  a  primitive  ambiguous  cycle  with  a  period  n.  Then, 


R  —  D(F0,  Fn) 

Proof:  Let  {F,:}  have  period  7 r  and  let  F0  and  Fn/ 2  be  the  two  symmetry  points  of  A. 
Then  F0*F0  =  1  =  Fn/2*  Fn/2,  with  Dpl  =  Dp2  —  0,  s  and  t  the  respective  first  coefficients. 
Therefore, 


0  —  D(F0  *  F0,  Fn/ 2  *  Fn/2)  —  2 D(F0,  Fn/2)  +  log [t/ s)  —  D(F0,  Fn) 

where  the  3rd  step  is  obtained  from  the  2nd  by  the  fact  that  the  product  in  D(F0,F, r/2) 
includes  the  last  denominator  t  and  not  the  first  denominator  s. 

Therefore,  D(Fq,F7T)  =  nR.  Considering  composition  of  Fq  with  forms  in  the  principal 
cycle,  clearly  D(F0,F^)  <  R,  so  that  D(F0,Fn)  =  R.  QED 

A. 6  Proof  of  Square  Forms  Factorization  (SQUFOF) 

It’s  not  certain  how  much  Shanks  may  have  rigorously  proven  concerning  distances,  but 
based  on  the  understanding  he  had  of  distance  and  infrastructure,  he  was  able  to  develop 
Square  Forms  Factorization.  A  short  example  will  demonstrate  and  explain  the  algorithm: 
let  N  =  3193.  Expanding  the  continued  fraction  (principal  cycle),  Qw  =  49.  The  quadratic 
form  for  this  is  F  —  49x2  +  58 xy  —  48 y2.  Since  49  is  a  perfect  square,  7x 2  +  58 xy  —  336 y2, 
which  reduces  with  Dp  =  0  to  G  =  7x2  +  100 xy  —  99 y2  is  a  quadratic  form  whose  square  is 
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F.  Therefore,  by  Theorem  A. 2. 8,  G  is  in  a  class  of  order  2  or  1,  so  that  G  is  an  ambiguous 
form,  so  that  there  are  two  points  of  symmetry  in  its  cycle.  Since  by  Theorem  A. 5. 2, 
2 D(GS,G)  =  D(1,F)  (mod  R).  So  D(GS,G)  =  D(l,F)/2  (mod  R/ 2).  Since  the  two  points 
of  symmetry  are  R/2  away  from  each  other,  this  means  that  there  is  a  symmetry  point  at 
distance  D(  1,  F)/2  behind  G.  Therefore,  a  point  of  symmetry  may  be  found  by  reversing  G 
and  traveling  this  short  distance.  Now  if  the  coefficient  at  this  symmetry  point  is  ±1,  then 
there  would  have  been  a  7  somewhere  before  F  in  the  continued  fraction  expansion.  If  the 
coefficient  is  2,  then  this  symmetry  point  could  be  composed  with  G  to  find  14  at  an  earlier 
point  in  the  principle  cycle.  Therefore,  the  symmetry  point  provides  a  nontrivial  factor  for 
N.  In  this  case,  after  6  steps  it  provides  31  as  a  factor  of  3193. 

The  second  phase  of  this  algorithm  can  be  made  significantly  (at  least  for  larger  numbers) 
faster  if  the  quadratic  forms  from  the  continued  fraction  expansion  with  indices  that  are 
powers  of  2  are  saved.  In  this  example,  F  =  Fw,  so  that  G  is  about  the  5th  form  in  its 
cycle11.  The  composition  of  G~l  with  F4  and  fq  is  close  and  a  simultaneous  search  in  both 
direction  from  there  quickly  finds  the  symmetry  point.  In  this  case,  it  is  only  necessary  to 
store  log2  k  forms  for  k  steps,  so  that  it  is  more  efficient  to  check  each  square  to  see  if  it 
works  than  to  check  each  square  root  against  the  previous  pseudo-squares  to  predict  whether 
it  will  work. 

Formally,  here  is  the  algorithm  for  factoring  N: 

11  Roughly,  since  in  this  case  5  «  6. 
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Qo  <-  1,  Po  <-  WN\,Qi  ^  N-P* 

r  <-  [v^J 

while  Q,  ^  perfect  square  for  some  %  even 


Pi  *-  biQi  -  Pi-i 

Qi+l  Qi- 1  +  bi(Pi- 1  —  Pi) 

if  i  =  2”  for  some  n 

_ _  p2  _  /y 

Store  (Qi,2  ■  Pt)  F0  =  (VQi,  2  •  P^,  ) 

Compose  Po  with  stored  forms  according  to  the 
binary  representation  of  i/2  and  store  result  to  Po- 

F0  =  (A,B,C) 

Qo  <-  \A\,P0<-B/2,Qi  <-  \C\ 

Qo  <-  Qi, Po  <-  PoPd  <-  Qo 
while  Pt  ^  Pi-i  and  ^  p^i 

Apply  same  recursive  formulas  to  (Qo,  Pq,Qi)  and  (qo,Po,Qi) 
If  Pi  =  Pi-i,  either  Qi  or  Qi/ 2  is  a  nontrivial  factor  of  N. 

If  pi  =  Pi-i,  either  qt  or  qi/2  is  a  nontrivial  factor  of  N. 


Finding  a  perfect  square  that  provides  a  factorization  is  the  slowest  part  of  the  algorithm, 
so  the  number  of  steps  required  to  obtain  this  is  a  good  measure  of  the  total  runtime.  Let  W 
be  the  number  of  forms  examined  before  a  square  for  is  found  that  provides  a  factorization. 
In  [23],  Shanks  states  that  for  N  having  k  distinct  prime  factors12, 


E(W)  =  ln(8) 


2  +  \[2  V N 
2  2k  —  2 


12Shanks  actually  let  N  have  k  +  1  distinct  prime  factors,  so  the  distance  looks  slightly  different  but  is 
equivalent. 
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Shanks  did  not  provide  a  proof,  and  it  actually  wasn’t  quite  right  for  all  N,  but  in  2004, 
Jason  Gower  analyzed  the  runtime  and  concluded: 

Conjecture  2.  [9]  Let  N  be  square- free  and  have  k  distinct  odd  prime  factors.  Let  W  be  the 
number  of  forms  that  SQLIFOF  must  examine  before  hireling  a  square  form  that  provides  a 
factorization.  Then, 


E(W)  ~ 


2(V2+l)tfNlog2  ijN  =  1  (mod  4^ 
ifN  =  2  or  3  (mod  4) 


(A. 37) 


Note  that  for  N  =  2  or  3  (mod  4),  this  is  equivalent  to  Shanks  estimate. 

Therefore,  Square  Forms  Factorization  probably  has  an  average  runtime  of  0(\/fN). 
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Appendix  B 
Source  Code 


B.l  Number  Theory  Functions 

This  first  program,  NumberTheory.c,  has  all  of  the  basic  functions  needed  for  composition 
of  quadratic  forms  and  cycling. 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  "gmp.h" 

struct  form  { 

mpz_t  QO ; 
mpz_t  P ; 
mpz_t  Q1 ; 
mpz_t  b; 

}; 

//prints  a  the  coefficients  of  a  quadratic  form  on  a  single  line 
void  outputF( struct  form  A) 

{ 

printf ("  "); 

mpz_out_str (stdout , 10 , A . QO) ; 
printf ("  "); 

mpz_out_str (stdout , 10 , A . P) ; 
printf ("  "); 
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mpz_out_str (stdout , 10,A.Q1) ; 
print f("  \n"); 

> 

//sets  dest  =  src 

void  ftof (struct  form  *dest, struct  form  src) 

{ 

mpz_set ( (*dest) . QO , src . QO) ; 
mpz_set ( (*dest) . P , src . P) ; 
mpz_set ( (*dest) . Q1 , src . Ql) ; 
mpz_set ( (*dest) .b, src .b) ; 

> 

//Calculates  modular  inverses.  Also  handles  if  A  and  base  are  not  relatively  prime. 
//Sets  d  =  gcd(A,base) .  Sets  inv  =  inverse  of  A  modulo  (base/d), 
void  inverse (mpz_t  inv,  mpz_t  d,  mpz_t  A,  mpz_t  base) 

{ 

mpz_t  t , one ; 
double  one_d  =  1; 
mpz_init (t) ; 

mpz_init_set_d(one , one_d) ; 
mpz_gcdext (d, inv,t , A, base) ; 
if (mpz_cmp(d,one) ) 

{ 

mpz_cdiv_q(A, A,d) ; 
mpz_gcdext (t , inv , t , A , base) ; 

> 

mpz_clear(t) ; 
mpz_clear (one) ; 

> 

//Solves  a  system  of  contruence  equations:  Chinese  remainder  theorem. 

//  Sets  soln  congruent  to  EQ[0]  modulo  EQ[1],  EQ  [2]  modulo  EQ [3] ,  etc. 
int  congruent (mpz_t  soln,  int  size,  mpz_t  EQ [] ) 

{ 

mpz_t  modu , cong , d , inv , temp , temp2 , prodSoFar ; 
int  i ; 

mpz_init (prodSoFar) ; 
mpz_init (modu) ; 
mpz_init (cong) ; 


mpz_init (d) ; 
mpz_init (inv) ; 
mpz_init (temp) ; 
mpz_init (temp2) ; 


mpz_set (soln,EQ [0] ) ; 
mpz_set (prodSoFar ,EQ [1] ) ; 
for(i  =  1;  i<size;  i++) 

{ 

mpz_set (cong,EQ [2*i] ) ; 
mpz_abs (modu,EQ  [2*i+l] ) ; 
inverse ( inv , d , prodSoF ar , modu) ; 
mpz_sub (temp , soln , cong) ; 
mpz_cdiv_r (temp2 ,  temp , d) ; 
if (mpz_sgn(temp2) ) 

{ 

mpz_clear (prodSoFar) ; 
mpz_clear (modu) ; 
mpz_clear (cong) ; 
mpz_clear (d) ; 
mpz_clear (inv) ; 
mpz_clear (temp) ; 
mpz_clear (temp2) ; 

return  0; 

} 

mpz_mul (temp , temp , inv) ; 
mpz_addmul (soln , temp , prodSoFar) ; 
mpz_cdiv_q(temp2 ,modu,d) ; 
mpz_mul (prodSoFar , prodSoFar , temp2) 
mpz_cdiv_r (soln, soln, prodSoFar) ; 

> 

mpz_clear (prodSoFar) ; 
mpz_clear (modu) ; 
mpz_clear (cong) ; 
mpz_clear(d) ; 
mpz_clear (inv) ; 
mpz_clear (temp) ; 
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mpz_clear (temp2) ; 

return  1 ; 

> 

//Takes  one  step.  Note  that  Qform  :=  [Q1 ,P1 ,Q2,bl] ,  all  positive, 
void  cFracStep (struct  form  *A,mpz_t  root) 

{ 

mpz_t  temp; 

mpz_init_set (temp, (*A) .P) ; 
mpz_add( (*A) .b,root , (*A) .P) ; 
mpz_f div_q( (*A) .b, (*A) .b, (*A) .Ql) ; 
mpz_submul ( (*A) .P, (*A) .b, (*A) .Ql) ; 
mpz_neg((*A) .P, (*A) .P) ; 
mpz_sub (temp, temp, (*A) .P) ; 
mpz_addmul ( (*A) .Q0,temp, (*A) .b) ; 
mpz_swap( (*A) .QO, (*A) .Ql) ; 

mpz_clear (temp) ; 

> 

//Tests  if  form  A  is  reduced. 

int  goodForm( struct  form  A,  mpz_t  root) 

{ 

mpz_t  IRoot ; 

mpz_init_set (IRoot ,root) ; 

mpz_sub (IRoot , root, A. Ql) ; 
mpz_abs (IRoot , IRoot) ; 

if ( (mpz_cmp(A.P,root)>0) I  I (mpz_cmp(A . P , IRoot) <0) ) 
return  0; 
return  1 ; 

> 

//Given  that  r  is  the  sqrt (square)  mod  base,  squares  base  and  finds  sqrt (square)  mod  new  base  congruent  to  r 
int  CompSquare (mpz_t  r,  mpz_t  square,  mpz_t  base) 

{ 

mpz_t  bs , inv,d, tempi ,temp2 ; 


mpz_init_set (bs ,base) ; 
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mpz_raul (base, base, base) ; 
if ( !mpz_cmp(base ,bs) ) 

{ 

mpz_clear (bs) ; 
return  1; 

> 

mpz_clear (bs) ; 

mpz_init (inv) ; 
mpz_init (d) ; 
mpz_init (tempi) ; 
mpz_init (temp2) ; 


inverse (inv, d,r ,base) ; 

mpz_f div_qr (tempi , temp2 , square , d) ; 

if (mpz_sgn (temp2) ) 

{ 

mpz_clear (inv) ; 
mpz_clear (d) ; 
mpz_clear (temp2) ; 
mpz_clear (tempi) ; 
return  0; 

> 


mpz_addmul(r, inv, tempi) ; 

if (mpz_tstbit (r , 0) ) 

mpz_sub(r ,r ,base) ; 

mpz_cdiv_q_ui(r ,r ,2) ; 
mpz_fdiv_r (r ,r ,base) ; 


mpz_clear (tempi) ; 
mpz_clear (temp2) ; 
mpz_clear (inv) ; 
mpz_clear(d) ; 
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return  1 ; 

> 

//Reverses  a  quadratic  form  and  takes  a  step  so  that  the  first  coefficient  is  unchanged, 
void  reverse (struct  form  *A,mpz_t  root) 

{ 

mpz_swap( (*A) . QO, (*A) . Ql) ; 
cFracStep(A,root) ; 

> 

//Reduces  A.  Does  not  use  the  original  value  of  A.Q1 
void  reduce (struct  form  *A,  mpz_t  root,  mpz_t  N) 

{ 

mpz_t  temp; 
mpz_init (temp) ; 

mpz_sub( (*A) .b,root , (*A) .P) ; 
mpz_f div_q( (*A) .b, (*A) .b, (*A) .QO) ; 
mpz_addmul ( (*A) .P, (*A) .b, (*A) .QO) ; 
mpz_mul(temp, (*A) .P, (*A) .P) ; 
mpz_sub( (*A) . Ql ,N,temp) ; 
mpz_cdiv_q( (*A) . Ql , (*A) . Ql , (*A) . QO) ; 

while ( !goodForm( (*A) ,root) ) 
cFracStep(A,root) ; 

mpz_clear (temp) ; 

> 

//Given  the  last  2  steps  in  either  the  denominator  or  the  numerater  of  the  continued  fraction  convergent 
//and  b,  calculates  the  next  value  and  shifts  R  forward  one  step, 
void  approx (mpz_t  R[2] ,  mpz_t  b) 

{ 

mpz_addmul (R [0] , b , R [1]  ) ; 
mpz_swap(R[0]  ,R[1]  )  ; 

> 

//Extended  gcd  algorithm.  Computes  d  =  gcd(x,y,z)  and  ax+by+cz  =  d. 
void  Xgcd3(mpz_t  d,mpz_t  a,  mpz_t  b,  mpz_t  c,  mpz_t  x,  mpz_t  y,  mpz_t  z) 
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{ 

mpz_t  d_xy , a_xy ; 
mpz_init (d_xy) ; 
mpz_init (a_xy) ; 

mpz_gcdext (d_xy ,  a ,  b ,  x ,  y) ; 
mpz_gcdext (d , a_xy , c , d_xy , z) ; 

mpz_mul (a, a, a_xy) ; 
mpz_mul (b , b , a_xy) ; 

mpz_clear (d_xy) ; 
mpz_clear (a_xy) ; 

> 

//Sets  res  =  A*B,  composition  of  quadratic  forms. 

void  compose (struct  form  *res,  struct  form  A,  struct  form  B,  mpz_t  root,mpz_t  N) 

{ 

mpz_t  tempi ,temp2,m,u,v,w; 

mpz_init (m) ; 

mpz_init (u) ; 

mpz_init (v) ; 

mpz_init (w) ; 

mpz_init (tempi) ; 

mpz_init (temp2) ; 

mpz_add(templ,A.P,B.P) ; 

Xgcd3(m,u,v,w,A.Q0,B.Q0,templ) ; 

mpz_mul (tempi , A . QO , B . P) ; 
mpz_mul (tempi , tempi ,u) ; 
mpz_mul (temp2 , A . P , B . QO) ; 
mpz_mul (temp2 , temp2 , v) ; 
mpz_add (tempi , tempi ,temp2) ; 
mpz_set (temp2,N) ; 
mpz_addmul (temp2 , A . P , B . P) ; 
mpz_mul(temp2,temp2 ,w) ; 
mpz_add (tempi , tempi ,temp2) ; 
mpz_cdiv_q( (*res) . P, tempi, m) ; 
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mpz_raul (tempi , A . QO , B . QO) ; 

mpz_mul (temp2 ,m,m) ; 

mpz_cdiv_q( (*res) .QO, tempi, temp2) ; 

reduce (res , root ,N) ; 


mpz_clear(m) ; 
mpz_clear(u) ; 
mpz_clear (v) ; 
mpz_clear(w) ; 
mpz_clear (tempi) ; 
mpz_clear (temp2) ; 

> 

B.2  Fast  Return  Functions 

FastReturn  was  the  name  Shanks  gave  to  the  variation  that  used  composition  to  test  a 
perfect  square  quickly.  This  program  provides  the  function  that  does  this. 

/*  FastReturn. h 

Used  by  SQUFOF  to  test  whether  a  number  is  square  and  use  saved  quadratic  forms  to  quickly  find  if  a  square 
//form  provides  a  factorization  of  N. 

*/ 

#include  "gmp.h" 

#include<math . h> 

int  subIsSquare(mpz_t  N,  mpz_t  root,  int  digit,  mpz_t  a) 

{ 

mpz_t  testN,newa,B; 
mpz_init_set (testN,N) ; 
mpz_submul (testN , a , a) ; 
if ( !mpz_sgn(testN) ) 

{ 

mpz_set (root , a) ; 
mpz_clear (testN) ; 


return  1 ; 
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> 

if (mpz_sgn(testN)<0) 

{ 

mpz_clear (testN) ; 
return  0; 

> 

if (mpz_tstbit (testN , 3*digit ) ) 

{ 

mpz_clear (testN) ; 
return  0; 

> 

int  b  =  mpz_tstbit (testN, 3*digit+l)+2*mpz_tstbit (testN, 3*digit+2) ; 
int  testa  =  mpz_tstbit (a,0)+2*mpz_tstbit (a, l)+4*mpz_tstbit (a, 2) ; 
b  =  (b*testa)#/«4; 

mpz_clear (testN) ;  // 

mpz_init_set_ui (B ,b) ; 
mpz_init (newa) ; 
mpz_mul_2exp(B,B,3*digit) ; 
mpz_add (newa ,a,B) ; 

int  t  =  subIsSquare(N, root ,digit+l , newa) ; 
if  (t) 

{ 

mpz_clear(B) ; 
mpz_clear (newa) ; 
return  1 ; 

> 

b  =  b+4; 

rapz_set_ui (B,b) ; 
mpz_mul_2exp(B,B,3*digit) ; 
mpz_add(newa, a,B) ; 

t  =  sublsSquare (N, root ,digit+l , newa) ; 
mpz_clear (B) ; 
mpz_clear (newa) ; 
return  t ; 


> 
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int  isSquare(mpz_t  N,mpz_t  root) 

{ 

mpz_t  IN; 
mpz_t  a; 
int  t; 

int  twos  =  0; 
double  i ; 

mpz_init_set (1N,N) ; 
while ( ! mpz_tstbit (IN , 0) ) 

{ 

twos++; 

mpz_fdiv_q_2exp(lN,lN, 1) ;  //Fix  this 

> 

if  (twos°/02) 
return  0; 
twos  =  twos/2; 
mpz_init (a) ; 

if (mpz_tstbit (IN, 1) I |mpz_tstbit (IN, 2) ) 
return  0; 

f or (i  =  1;  i<8 ;  i=i+2) 

mpz_set_d(a, i) ; 
t  =  subIsSquare(lN,root,l,a) ; 
if  (t) 

{ 

mpz_clear(a) ; 
mpz_clear (IN) ; 

rapz_mul_2exp (root , root , twos) ; 
return  1; 

> 

> 

mpz_clear (a) ; 
mpz_clear (IN) ; 
return  0; 


} 
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//Given  that  the  first  coefficient  of  A  has  sqrt  =  r  and  index  i.  saved[]  is  a  list  of  previous  forms 
//such  that  saved[i]  =  F_{30*2~i} .  root  =  sqrt(N),  where  N  is  the  discriminant  (the  number  being  factored), 
//tr  is  set  to  1  if  a  factor  is  found. 

int  testF(struct  form  A, struct  form  saved [],int  i,mpz_t  r,mpz_t  root,mpz_t  N,int  *tr) 

i 

struct  form  B,myA; 
mpz_t  BP,AP,d; 
i  =  i/2; 
int  size  =  24; 
int  os; 

int  pow  =  503316480;  //30*2~24 

mpz_init_set(myA.Q0,A.Q0) ;  mpz_init_set (myA . P , A . P) ;  mpz_init_set(myA.Ql,A.Ql) ;  mpz_init (myA .b) ; 

reverse (&myA, root) ; 

mpz_set (myA . Q0 , r ) ; 

reduce (&myA, root ,N) ; 

while  (size  >=  0) 

{ 

if (i>=pow) 

{ 

i  =  i-pow; 

compose (&myA, myA, saved [size] ,root,N) ; 

> 

size  =  size-1; 
pow  =  pow/2; 

> 

mpz_init (B . Q0) ;  mpz_init (B . P) ;  mpz_init (B . Ql) ;  mpz_init (B .b) ; 
ftof (&B ,myA) ; 
reverse (&B , root) ; 

mpz_init (d) ; 
mpz_init_set_ui (BP , 1) ; 
mpz_init_set_ui (AP , 1) ; 

if ( !mpz_cmp(myA.P,B.P)) 

{ 

mpz_gcd(d,myA.P,N) ; 

mpz_clear (BP) ;  mpz_clear (AP) ; 


mpz_clear (B . QO) ;  mpz_clear (B .P) ;  mpz_clear (B . Ql) ;  mpz_clear (B .b) ; 
mpz_clear (myA . QO) ;  mpz_clear(myA.P) ;  mpz_clear (myA. Ql) ;  mpz_clear(myA.b) 

if (mpz_cmp_ui (d ,  2)  >0) 

{ 

mpz_out_str (stdout , 10 ,d) ; 
mpz_clear(d) ; 

(*tr)  =  1; 
return  1; 

> 

else 

{ 

mpz_clear(d) ; 

(*tr)  =  0; 
return  0; 

> 

> 

os  =  0; 

while  (mpz_cmp (myA . P , AP) &&mpz_cmp (B . P , BP) ) 

{ 

mpz_set (AP ,myA . P) ; 
mpz_set (BP , B . P) ; 
cFracStep(&myA,root) ; 
cFracStep(&B,root) ; 
os  =  os+1; 

> 

if  ( !mpz_cmp(myA.P,AP)) 

mpz_gcd(d,myA. Q0,N) ; 

mpz_clear (BP) ;  mpz_clear (AP) ; 

mpz_clear (B . Q0) ;  mpz_clear (B .P) ;  mpz_clear (B . Ql) ;  mpz_clear (B .b) ; 
mpz_clear (myA . Q0) ;  mpz_clear(myA.P) ;  mpz_clear (myA. Ql) ;  mpz_clear (myA . b) 

if (mpz_cmp_ui (d , 2) >0) 

{ 


mpz_out_str (stdout , 10, d) ; 
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printf ("\n") ; 
mpz_clear (d) ; 
(*tr)=l; 
return  1; 

> 

else 

{ 

mpz_clear(d) ; 
(*tr)=0; 
return  0; 

> 


> 

if ( ! mpz_cmp (B . P , BP) ) 

{ 

mpz_gcd(d,B . Q0 ,N) ; 

mpz_clear (BP) ;  mpz_clear(AP) ; 

mpz_clear (B . Q0) ;  mpz_clear (B .P) ;  mpz_clear (B . Ql) ;  mpz_clear (B .b) ; 
mpz_clear(myA. Q0) ;  mpz_clear(myA.P) ;  mpz_clear (myA. Ql) ;  mpz_clear(myA.b) ; 

if (mpz_cmp_ui (d , 2) >0) 

{ 

mpz_out_str (stdout , 10 ,d) ; 
printf ("\n") ; 
mpz_clear(d) ; 

(*tr)=l ; 
return  1; 

> 

else 

{ 

mpz_clear(d) ; 

(*tr)=0; 
return  0; 

> 

> 

mpz_clear (BP) ;  mpz_clear(AP) ; 

rapz_clear (B . Q0) ;  mpz_clear (B .P) ;  mpz_clear (B . Ql) ;  mpz_clear (B .b) ; 


mpz_clear (myA. QO) ;  mpz_clear(myA.P) ;  mpz_clear (myA. Ql) ;  mpz_clear (myA.b) ; 
return  0; 


no 


> 


B.3  Segments  Parallelization 

This  code,  parSQUFOFsegs.c,  uses  segments  to  distribute  the  cycle  of  quadratic  forms  be¬ 
tween  multiple  processors.  It  is  designed  for  up  to  64  processors,  although  it  would  be  easy 
to  modify  for  a  larger  system  or  for  distributed  processing.  It  reads  from  standard  output. 
Each  processor  also  records  some  basic  information  about  the  factors  it  finds  to  a  hie  named 
“lSXX.txt”,  where  “XX”  is  the  two  digit  process  number.  At  the  end  of  the  program,  the 
master  compiles  all  output  into  the  hie  “SegmentsData.txt”.  It  has  several  other  options. 

-n=x  specifies  that  there  will  be  x  integers  to  factor,  -d  specifies  that  some  basic  informa¬ 
tion  should  be  output  for  each  factor  found.  -D  specifies  that  thorough  information  should 
be  provided  about  the  progress  of  each  processor  to  assist  in  debugging,  -i  specifies  that  the 
computer  should  be  interactive  with  the  user 

Thus,  a  standard  format  for  the  command  would  be: 
mpirnn  -np  50  parSQUFOFsegs  jnumbers.txt  -n=100  -d 

which  would  use  50  processors  to  factor  100  integers  from  the  hie  numbers.txt  and  output 
basic  information  about  each  to  the  screen. 

//parSQUFOFsegs . c 
//mm  1/C  Steve  McMath 
//April,  2005 


#include  "mpi.h" 

#include  <stdio.h> 

#include  "NumberTheory .h" 

#include  "/usr/local/include/gmp .h" 
#include  <time.h> 


#include  "FastReturn.h1 


Ill 


#include  "gmpmpi_segments .h" 

int  main (int  argc,  char  *argv[]) 

{ 

const  int  STOP  =  10; 

const  int  finished  =  11; 

const  int  segSize  =  25; 

const  int  Npassl  =  12; 

const  int  Npass2  =  13; 

char  f name [9]  =  "lS00.txt"; 

int  myrank,i, j ,Ni ,np , flag , done , tens , single ; 

mpz_t  N, root ,firstBack, factor ; 

struct  form  Start; 

struct  form  back  [segSize] ; 

char  *N_str; 

int  lengthN; 

int  dDi [3]  =  {0,0,0}; 

int  maxNi  =  100; 

MPI_Status  status; 
time_t  start, endt; 
double  dif ; 

FILE  *0UT, *last0UT, *IN ; 

MPI_Init (&argc ,  &argv) ; 

MPI_Comm_rank (MPI_C0MM_W0RLD ,  &myrank) ; 
MPI_Comm_size (MPI_C0MM_W0RLD ,  &np) ; 

//Idiot  check 

if  ((np  <  2)  II  (np  >  50)) 

{ 

MPI_Finalize() ; 
return  0; 

> 

//Read  in  parameters  from  user. 
for(i  =  1;  i<argc;  i++) 

{ 

if (argv[i]  [0]  ==  ’-’) 

{ 
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if (argv[i] [1]  ==  ’n’) 

{ 

maxNi  =  0; 

for(j  =  3;  argv[i]  [j]  !=NULL;  j++) 

{ 

maxNi  =  10*maxNi+(argv [i] [j] - ’O' ) ; 

} 

} 

if (argv[i] [1]  ==  ’d’) 
dDi  [0]  =  1; 
if (argv[i] [1]  ==  ’D’) 

{ 

dDi  [0]  =  1; 
dDi  [1]  =  1; 

} 

if  (argv[i]  [1]  ==  ’i’) 
dDi  [2]  =  1; 

> 

> 

//Set  2  digits  in  filename  to  processor  number. 
fname[2]  =  ’O’+myrank/lO; 
f name  [3]  =  ’  0’+myrank"/,10; 

OUT  =  f open(f name , "w" ) ; 

for(Ni=0;  Ni  <  maxNi;  Ni++) 

{ 


f flush (OUT) ; 
flag  =  0; 
done  =  0; 
mpz_init (N) ; 

if (myrank==0)  //Read  in  N  and  pass  it  to  the  other  processors. 

{ 

fprintf  (OUT,  "Starting  */,i\n"  ,Ni)  ; 
mpz_inp_str (N, stdin, 10) ; 
lengthN  =  l+mpz_sizeinbase(N,36) ; 

N_str  =  malloc(lengthN*sizeof (char) ) ; 
mpz_get_str(N_str ,36,N) ; 
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for(i  =  1;  i<np;  i++) 

{ 

MPI_Send(&lengthN, 1 ,MPI_INT, i ,Npassl ,MPI_COMM_WORLD) ; 

MPI_Send(N_str , lengthN , MPI_CHAR , i , Npass2 , MPI_COMM_WORLD) ; 

> 

free(N_str) ; 
time(&start) ; 

> 

else  //Receive  N  from  the  master. 

{ 

MPI_Recv(&lengthN , 1 , MPI_INT , 0 , Npass 1 , MPI_C0MM_W0RLD , &status) ; 

N_str  =  malloc(lengthN*sizeof (char) ) ; 

MPI_Recv(N_str , lengthN , MPI_CHAR , 0 , Npass2 , MPI_C0MM_W0RLD , &status) ; 
mpz_set_str (N,N_str ,36) ; 
free(N_str) ; 
time(&start) ; 

> 

//Initialize  memory  and  calculate  first  form. 
mpz_init (root) ; 
mpz_sqrt (root ,N) ; 

mpz_init_set_ui (Start . QO, 1) ; 
mpz_init_set (Start .P,root) ; 
mpz_init_set (Start . Q1 ,N) ; 
mpz_submul (Start . Ql, root, root) ; 
mpz_init (Start . b) ; 

//Take  30  steps  forward  from  the  start. 
mpz_init (f irstBack) ; 
for(i  =  0;  i<30;  i++) 

{ 

cFracStep(& (Start) ,root) ; 
cFracStep(&(Start) ,root) ; 

if (isSquare (Start . QO , f irstBack) &&testF (Start , back , Start , 0 , f irstBack , root , N , &done , factor , dDi) ) 

{ 

mpz_clear (root) ; 
mpz_clear (Start . QO) ; 
mpz_clear (Start .P) ; 
mpz_clear (Start . Ql) ; 
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mpz_clear (Start .b) ; 
mpz_clear(N) ; 

MPI_Finalize() ; 
return  0; 

> 

> 

//Compute  forms  used  by  Fast  Return, 
ford  =  0;  i<segSize;  i++) 

{ 

mpz_init_set (back [i] . Q0 , Start . Q0) ; 
mpz_init_set (back [i] . P , Start . P) ; 
mpz_init_set (back [i] .  Ql , Start . Q 1 ) ; 
mpz_init (back [i] .b) ; 
compose (&Start , Start , Start , root ,N) ; 

> 

if  (myrank  ==0) 

{ 

int  found [64]; 

MPI_Request  waiting [64] , stop [64] ; 

struct  segment  Next ; 

struct  form  hlncrement ,nextHstart ; 

ford  =  1;  i<np;  i++)  //Receive  messages  from  ready  processors. 

{ 

waiting [i]  =  0; 

MPI_Irecv (&f ound [i] , 1 , MPI.INT , i finished , MPI_C0MM_W0RLD , &wait ing [i] ) ; 

> 

//hlncrement  is  used  as  the  first  step  before  Fast  Return. 
mpz_init_set (hlncrement .Q0, back [segSize-1] .Q0) ; 
mpz_init_set (hlncrement . P ,back [segSize-1] .P) ; 
mpz_init_set (hlncrement . Q1 , back [segSize-1]  . Ql) ; 
mpz_init (hlncrement .b) ; 

//Next. start  is  the  first  form  in  the  segment. 
mpz_init_set (Next. start .Q0, back [0] .Q0) ; 
mpz_init_set (Next . start . P , back [0] . P) ; 
mpz_init_set (Next . start . Ql ,back[0] . Ql) ; 
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mpz_init (Next . start . b) ; 

//Next. end  is  the  last  form  in  the  assigned  segment. 
mpz_init_set (Next . end . QO , Start . QO) ; 
mpz_init_set (Next . end . P , Start . P) ; 
mpz_init_set (Next . end. Q1 , Start . Ql) ; 
mpz_init (Next . end . b) ; 

//half Start  is  used  in  the  first  step  before  Fast  Return. 
mpz_init_set (Next . half Start . QO , back [0] . QO) ; 
mpz_init_set (Next .half Start .P, back [0] .P) ; 
mpz_init_set (Next .half Start . Ql ,back [0] . Ql) ; 
mpz_init (Next . half Start . b) ; 

mpz_init_set (nextHstart . QO , hlncrement . QO) ; 
mpz_init_set (nextHstart . P , hlncrement . P) ; 
mpz_init_set (nextHstart . Ql , hlncrement . Ql) ; 
mpz_init (nextHstart .b) ; 

//Records  which  segment  this  is. 

Next.segnum  =  0; 
i  =  1; 

while(! done)  //Checks  for  free  processors,  sends  them  a  segment,  then  calculates  the  next. 

{ 

if (i>=np) 
i  =  1; 
flag  =  0; 

MPI_Test (&waiting [i] ,&f lag,&status) ; 
if (flag) 

{ 

if (found [i] ) 

{ 

if (dDi  [1] ) 

printf ("Something  found  by  rank  °/0d\n",i); 
done  =  found [i] ; 

> 

else 

{ 

Send(Next , i) ; 
if (dDi  [1]) 
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printf  ("Segment  °/0d  given  to  rank  0/0d\n"  ,Next .  segnum,  i)  ; 

Next . segnum++; 

MPI_Irecv(&f ound [i] , 1 , MPI_INT , i finished , MPI_C0MM_W0RLD , &waiting [i] ) ; 

f tof (&Next . start , Next . end) ; 

ftof (&Next .half Start ,nextHstart) ; 

compose (&nextHstart , nextHstart , hlncrement , root ,  N)  ; 

compose (&Next . end , nextHstart , nextHstart , root , N) ; 

> 

> 

i++ ; 

} 

//Let  all  the  processors  know  when  one  finds  a  factor. 
for(i=l;  i<np;  i++) 

MPI_Isend (&done , 1 , MPI.INT , i , STOP , MPI_C0MM_W0RLD , &stop [i] ) ; 

//Clear  memory. 

mpz_clear (hlncrement . QO) ;  mpz_clear (hlncrement .P) ;  mpz_clear (hlncrement . Ql) ;  mpz_clear (hlncrement .b) ; 
mpz_clear (Next . start . QO) ;  mpz_clear (Next . start .P) ;  mpz_clear (Next . start . Ql) ;  mpz_clear (Next . start .b) ; 
mpz_clear (Next . end. QO) ;  mpz_clear (Next . end. P) ;  mpz_clear (Next .end. Ql) ;  mpz_clear (Next . end. b) ; 
mpz_clear (Next .half Start . QO) ;  mpz_clear (Next .half Start .P) ;  mpz_clear (Next .half Start . Ql) ; 
mpz_clear (Next . half Start . b) ; 

mpz_clear (nextHstart . QO) ;  mpz_clear (nextHstart .P) ;  mpz_clear (nextHstart . Ql) ;  mpz_clear (nextHstart .b) ; 

//Wait  for  all  processors  to  get  the  message. 
for(i=l;  i<np;  i++) 

MPI_Wait (&waiting [i] ,&status) ; 

> 

else 

{ 

struct  segment  Mine; 

MPI_Request  last; 

MPI_Request  ready; 

long  long  int  index; 

int  test=0, found  =  0,end  =  0; 

//Initialize  memory. 

mpz_init (Mine . start . QO) ;  mpz_init (Mine . start . P) ;  mpz_init (Mine . start . Ql) ;  mpz_init (Mine . start .b) ; 
mpz_init (Mine . end. QO) ;  mpz_init (Mine. end. P) ;  mpz_init (Mine . end. Ql) ;  mpz_init (Mine . end. b) ; 
mpz_init (Mine .half Start . QO) ;  mpz_init (Mine .half Start .P) ;  mpz_init (Mine .half Start . Ql) ; 
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mpz_init (Mine . half Start . b) ; 

//Set  up  to  receive  the  message  when  a  processor  finishes. 

MPI_Irecv (&done , 1 , MPI.INT , 0 , STOP , MPI_C0MM_W0RLD , ftlast ) ; 
if (dDi [1] ) 

printf  ("°/0d,  done  =  °/0d  recieved\n"  ,myrank,done)  ; 
while ((! flag) &&(! end) )  //Each  cycle  through  this  loop  receives  a  segment  from  master. 

{ 

//Send  ’ready’  to  master. 

MPI_Isend (&f ound , 1 , MPI_INT , 0  finished , MPI_C0MM_W0RLD , &ready) ; 

//Receive  a  segment. 

Recv(&Mine) ; 

if ( (Mine . segnum==(np-2) )&&dDi [1] ) 

i 

time(&endt) ; 

dif  =  difftime(endt , start) ; 

printf  ("All  processes  started  in  °/0e  seconds .  \n"  , dif ) ; 

> 

reduce (&Mine . start , root ,  N)  ; 

index  =  0; 

//Basic  step.  Checks  for  messages  every  cycle. 

//Does  Fast  Return  if  it  finds  a  square. 

while ( ( ( ! isSquare (Mine . start . Q0 ,f irstBack) ) I  I ( ! testF (Mine . start , back, Mine .half Start , index, firstBack, 
root ,N,&test , factor , dDi) ))&&(! end) ) 

{ 

cFracStep (&Mine . start , root) ; 
cFracStep(&Mine . start , root) ; 

//Check  for  Start 

if ( ( !mpz_cmp_ui (Mine . start . Q0 , 1) ) I  I ( !mpz_cmp_ui (Mine . start . Q1 , 1) ) ) 

{ 

end  =  1; 
found  =  -1; 
if (dDi  [0] ) 

printf  ("Start  found  in  segment  °/0d"  ,Mine  .  segnum)  ; 

> 


//Check  for  end  of  segment. 

if ( ( ! mpz_cmp (Mine . start . QO , Mine . end . QO) ) && ( ! mpz_cmp (Mine . start . P , Mine . end . P) ) ) 
end  =  1 ; 

//Check  for  message. 

MPI_Test (&last ,&f lag,&status) ; 
if (flag) 

{ 

if (done==l) 
end  =  1; 

> 

index++; 

> 

if  (test) 

{ 

found  =  1; 
end  =  1; 
time(&endt) ; 

dif  =  difftime(endt , start) ; 
if (dDi  [0] ) 

printf  ("#/«i ,  Rank  #/*i,  0/0e,  #/«i,  " ,Ni ,myrank, dif , Mine . segnum) ; 
fprintf  (OUT ,  "#/«i ,  Rank  °/#i,  °/0e,  °/ei,  " ,Ni ,myrank, dif , Mine . segnum) ; 
mpz_out_str (stdout , 10, factor) ; 
printf ("\n") ; 
fprintf (OUT, "\n") ; 

> 

MPI_Test (&last ,&f lag,&status) ; 

> 

MPI_Isend (&f ound , 1 , MPI_INT , 0  finished , MPI_C0MM_W0RLD , &ready) ; 

mpz_clear (Mine . start . QO) ;  mpz_clear (Mine . start .P) ;  mpz_clear (Mine . start . Ql) ;  mpz_clear (Mine . start .b) 
mpz_clear (Mine . end. QO) ;  mpz_clear (Mine . end. P) ;  mpz_clear (Mine .end. Ql) ;  mpz_clear (Mine . end. b) ; 
mpz_clear (Mine .half Start . QO) ;  mpz_clear (Mine .half Start .P) ;  mpz_clear (Mine .half Start . Ql) ; 
mpz_clear (Mine . half Start . b) ; 

> 

time(&endt) ; 

dif  =  dif f time (endt , start) ; 
if (dDi [1]) 

printf  ("Rank  =  °/0d:  Done  in  °/0e  sec\n"  ,myrank, dif ) ; 
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//Clear  memory  used. 
for(i=0;  i<segSize;  i++) 

{ 

mpz_clear (back[i]  .QO) ; 
mpz_clear (back[i] .P) ; 
mpz_clear (back[i]  .Ql) ; 
mpz_clear (back[i] .b) ; 

> 

mpz_clear (Start . QO) ; 
mpz_clear (Start . P) ; 
mpz_clear (Start . Ql) ; 
mpz_clear (Start .b) ; 
mpz_clear (N) ; 
mpz_clear (root) ; 

> 

f close (OUT) ; 

if (myrank==0)  //Compile  the  information  from  the  different  processors  into  a  single  file. 

{ 

lastOUT  =  f open("SegmentsData. txt" , "w") ; 

for(i=0;  i<np;  i++) 

{ 

tens  =  i/10; 
f name [2]  =  ’ O’+tens; 
f name  [3]  =  ’  0’+i#/010; 

IN  =  f open(fname , "r") ; 
single  =  fgetc(IN); 
do 
{ 

fprintf  (lastOUT,  "#/0c"  , single)  ; 
single  =  fgetc(IN); 

}  while (single !=E0F) ; 
fprintf (lastOUT, "\n") ; 
f close(IN) ; 
f flush (lastOUT) ; 


> 

f close (lastOUT) ; 
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> 


MPI_Finalize() ; 
return  0; 


> 


B.3.1  Segments  Communication  Functions 

This  file,  gmpmpi_segments.li,  is  used  by  parSQUFOFsegs.c  to  send  the  data  related  to  a 
segment  from  the  master  processor  to  the  other  processors. 


struct  segment  { 

struct  form  start; 
struct  form  end; 
struct  form  half Start; 
int  segnum; 


void  Send(struct  segment  A,  int  dest) ; 
void  Recv(struct  segment  *A)  ; 


//Determines  the  lengths  of  each  intiger  in  base  36  and  allocates  space  for  them  to  be  stored  as  integers.  The 
//lengths  can  be  sent  to  the  other  processor  to  allocate  space  and  be  able  to  receive, 
void  form_type (struct  segment  A,  char  *A_str[6],  int  lengths [6] ) 

{ 

int  i ; 


lengths [0] 
lengths [1] 
lengths [2] 
lengths [3] 
lengths [4] 
lengths [5] 


l+mpz_sizeinbase (A . start . Q0 , 36) ; 
l+mpz_sizeinbase (A . start . P , 36) ; 
l+mpz_sizeinbase (A . end . Q0 , 36) ; 
l+mpz_sizeinbase (A . end . P , 36) ; 
l+mpz_sizeinbase (A . half Start . Q0 , 36) ; 
l+mpz_sizeinbase (A . half Start . P , 36) ; 


ford  =  0;  i<6;  i=  i+1) 

{ 


A_str[i]  =  malloc (lengths [i] *sizeof (char) ) ; 
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> 

mpz_get_str (A_str [0] , 36 , A . start . Q0) ; 
mpz_get_str (A_str [1] , 36 , A . start . P) ; 
mpz_get_str (A_str [2] , 36 , A . end . Q0) ; 
mpz_get_str (A_str [3] , 36 , A . end . P) ; 
mpz_get_str (A_str [4] ,36, A. half St art .Q0) ; 
rapz_get_str (A_str [5] , 36 , A . half Start . P) ; 

> 

//allocates  space  for  the  lengths  sent 

void  ref orm_type( char  *A_str [6] , int  lengths [6] ) 

{ 

int  i ; 

for(i  =  0;  i<6;  i++) 

A_str[i]  =  malloc (lengths [i] *sizeof (char) ) ; 

> 

//After  receiving  the  strings,  converts  them  to  mpz  and  stores  them  correctly, 
void  unpack (struct  segment  *A,  char  *A_str[6]) 

{ 

mpz_set_str ( (*A) . start . Q0 , A_str [0] , 36) ; 
mpz_set_str ( (*A) . start .P,A_str [1] ,36) ; 
mpz_set_str ( (*A) . end . Q0 , A_str  [2] , 36) ; 
mpz_set_str ( (*A) . end.P, A_str  [3] ,36) ; 
mpz_set_str ( (*A) . half Start . Q0 , A_str [4] , 36) ; 
mpz_set_str ( (*A) .half St art .P,A_str [5] ,36) ; 

> 

//Sends  a  segment  A  to  processor  #dest 
void  Send (struct  segment  A,  int  dest) 

{ 

char  *As  [6] ; 
int  lengths [6] ; 

const  int  tags [6]  =  {1,2, 3, 4, 5, 6}; 
const  int  extra  =  7 ; 
int  i; 


f orm_type (A , As , lengths) ; 
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MPI_Send(&lengths [0] , 6 , MPI_INT , dest , 0 ,MPI_CQMM_W0RLD) ; 
for(i  =  0;  i<6;  i++) 

MPI_Send(As [i] , lengths [i] ,MPI_CHAR, dest , tags [i] ,MPI_C0MM_W0RLD) ; 

MPI.Send (& (A . segnum) , 1 , MPI_INT , dest , extra , MPI_C0MM_W0RLD) ; 
for(i  =  0;  i<6;  i++) 
free(As [i] ) ; 

> 

//Receive  segment  A  from  processor  0 
void  Recv (struct  segment  *A) 

{ 

char  *As [6] ; 
int  lengths [6] ; 

const  int  tags [6]  =  {1,2, 3, 4, 5, 6}; 
int  i; 

MPI_Status  halt; 
const  int  extra  =  7 ; 

MPI.Recv (&lengths [0] , 6 , MPI.INT , 0 , 0 , MPI_C0MM_W0RLD , fehalt) ; 
ref orm_type (As , lengths) ; 
ford  =  0;  i<6;  i++) 

MPI_Recv(As [i] , lengths [i] ,MPI_CHAR,0,tags [i] ,MPI_C0MM_W0RLD,&halt) ; 
unpack (A, As) ; 

MPI_Recv(& ( (*A) . segnum) , 1 , MPI.INT , 0 , extra , MPI_C0MM_W0RLD , fehalt) ; 
for(i=0;  i<6;  i++) 
free(As  [i] ) ; 


B.4  Multipliers  Parallelization 

This  code,  parSQUFOFmults.c,  uses  different  multipliers  to  distribute  the  factorization  al¬ 
gorithm  between  multiple  processors.  It  is  currently  configured  for  up  to  50  processors, 
although  it  would  be  easy  to  modify  it  for  a  larger  system  or  for  distributed  processing. 
Commands,  input,  and  output  are  the  same  as  parSQUFOFsegs,  except  that  the  hies  are 
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named  “2SXX.txt”  and  the  final  file  output  is  named  “MultipliersData.txt”. 


//parSQUFOFmults . c 
//mm  1/C  Stephen  McMath 
//May  2005 


#include  "mpi.h" 

#include  <stdio.h> 

#include  "NumberTheory .h" 

#include  "/usr/local/include/gmp .h" 

#include  <time.h> 

#include  "FastReturn.h" 

#define  NUMPROC  51  //Number  of  processors  configured  for. 

//Single  implementation  of  SQUFOF. 

int  SQUF0F(mpz_t  N,mpz_t  factor, int  dDi [3] ) ; 

int  main (int  argc,  char  *argv[]) 

{ 

//Constants  for  sending  message  to  stop  working. 

const  int  ST0P[64]  =  {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31, 
32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63}; 
const  int  Npassl  =  12;  //Constants  for  passing  N  from  master  to  slaves, 
const  int  Npass2  =  13; 
int  myrank,np,i,j , success; 

int  dDi [3]  =  {0,0,0};  //d  for  basic  information.  D  for  thorough  debugging.  i  for  interactive 
int  maxNi  =  100;  //Can  also  be  changed  by  user. 
mpz_t  N,  factor,  temp; 

char  *N_str;  //Used  for  passing  N  from  the  master, 
int  lengthN; 

char  fname[9]  =  "2S00.txt";  //Stores  basic  information  from  each  processor. 

MPI_Request  stop [NUMPROC] ; 

//Multiples  chosen  to  minimize  size  of  largest  prime  used  and  arranged  in  increasing  order, 
int  mults [NUMPROC]  =  {1,2,3,5,6,7,10,11,13,14,15,21,22,26,30,33,35,39,42,55,65,66,70,77,78,91,105,110,130, 
143,154,165,182,195,210,231,273,286, 330 , 338 , 385 , 390 , 429 , 455 ,462,546,715, 770 ,858,1155, 2310} ; 

MPI_Status  status; 
time_t  start, endt; 
double  dif ; 


int  Ni,tens; 
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int  single; 

FILE  *OUT; 

FILE  *IN,  *lastOUT ; 

time(festart) ; 

MPI_Init (&argc ,  &argv) ; 

MPI_Comm_rank(MPI_COMM_ WORLD,  &myrank) ; 

MPI_Comm_size (MPI_C0MM_W0RLD,  Stop) ; 

//Idiot  check 

if ( (np  <  2)  ||  (np  >  NUMPROC) ) 

{ 

MPI_Finalize () ; 
return  0; 

> 

//Read  in  parameters  from  user, 
ford  =  1;  i<argc;  i++) 

{ 

if  (argv[i]  [0]  == 

{ 

if (argv[i] [1]  ==  ’n’) 

{ 

maxNi  =  0; 

for(j  =  3;  argv [i]  [j]  !  =NULL;  j++) 

{ 

maxNi  =  10*maxNi+(argv [i] [j] - ’O’ ) ; 

} 

} 

if (argv [i]  [1]  ==  ’d’) 
dDi  [0]  =  1; 
if (argv [i]  [1]  ==  ’D’) 

{ 

dDi  [0]  =  1; 
dDi  [1]  =  1; 

} 

if (argv [i]  [1]  ==  ’i’) 
dDi  [2]  =  1; 

> 
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> 

mpz_init (N) ; 

//Set  2  digits  in  filename  to  processor  number, 
tens  =  myrank/10; 
f name [2]  =  ’O’+tens; 
f name  [3]  =  ’O’+myrank/olO; 

OUT  =  f open (f name, "w") ; 

for(Ni=0;  Ni<maxNi;  Ni++) 

{ 


f flush (OUT) ; 

if (myrank==0)  //Read  in  N  and  pass  it  to  the  other  processors. 

{ 

if (dDi [2] ==1) 
printfO'N:  \n"); 
mpz_inp_str (N , stdin , 10) ; 
lengthN  =  l+mpz_sizeinbase(N,36) ; 

N_str  =  malloc (lengthN*sizeof (char) ) ; 
mpz_get_str (N_str ,36,N) ; 
for(i  =  1;  i<np;  i++) 

{ 

MPI_Send(&lengthN, 1 ,MPI_INT, i ,Npassl ,MPI_C0MM_W0RLD) ; 

MPI_Send (N_str , lengthN , MPI.CHAR , i , Npass2 , MPI_C0MM_W0RLD) ; 

> 

free(N_str) ; 

fprintf  (OUT,  "Starting  Ni  =  °/0d\n",Ni); 
time(&start) ; 

> 

else  //Receive  N  from  the  master. 

{ 

MPI_Recv(&lengthN , 1 , MPI_INT , 0 , Npass 1 , MPI_C0MM_W0RLD , &status) ; 
N_str  =  malloc (lengthN*sizeof (char) ) ; 

MPI_Recv(N_str , lengthN , MPI_CHAR , 0 , Npass2 , MPI_C0MM_W0RLD , &status) ; 
mpz_set_str (N,N_str ,36) ; 
free(N_str) ; 
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time(&start) ; 

> 

//Each  processor  multiplies  N  by  its  multiple. 
mpz_mul_ui (N,N,mults [myrank] ) ; 

if  (dDi[l]) 

print  f  ("Rank  °/0d,  *°/0d\n"  ,  my  rank,  mult  s  [myrank]  )  ; 

mpz_init (factor) ; 

//Apply  standard  SQUFOF  algorithm  (with  Fast  Return)  to  the  new  N. 
success  =  SQUFOF (N,f actor ,dDi) ; 

if (success)  //Output  information  according  to  user  parameters. 

{ 

mpz_init (temp) ; 

mpz_gcd_ui(temp,f actor,mults [myrank] ) ; 
mpz_fdiv_q(f actor , factor , temp) ; 
mpz_clear (temp) ; 
if (dDi [1] ) 

printf  ("Rank  °/0d,  success  !  \n"  , myrank)  ; 
time(&endt) ; 

dif  =  difftime(endt , start) ; 

f  printf  (OUT,  "°/0d,  °/0d,  °/0e,  "  ,Ni ,  mult  s  [myrank]  ,  dif  )  ; 
if (dDi  [0] ) 

{ 

printf  ("Ni:  °/0d,  mult:  °/0d,  time:  °/0e  sec,  factor:  "  ,Ni , mult s [myrank]  , dif )  ; 
mpz_out_str (stdout , 10, factor) ; 
printf ("\n") ; 

} 

mpz_out_str (OUT, 10, factor) ; 
f printf (OUT, "\n") ; 

for(i  =  0;  i<np;  i++)  //Inform  the  other  processors. 

{ 

if (i ! =myrank) 

MPI_Isend(&success,l,MPI_INT,i, STOP [myrank] ,MPI_C0MM_W0RLD,&stop [i] ) ; 

> 

for(i=0;  i<np;  i++)  //Wait  for  each  processor  to  receive  the  message. 
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{ 

if (i ! =myrank) 

MPI_Wait (&stop [i] ,&status) ; 

> 

> 

time(&endt) ; 
if  (dDi[l]) 

{ 

dif  =  difftime(endt , start) ; 

printf("Rank  °/0d,  done  in  70e\n"  ,myrank,dif )  ; 

> 

if (myrank ! =0)  //Inform  master  when  ready  for  next  N. 

MPI_Isend(&success , 1 ,MPI_INT,0,-1 ,MPI_C0MM_W0RLD,&stop [0] ) ; 
else 

for(i=l;  i<np;  i++) 

MPI_Recv(&success , 1 ,MPI_INT, i ,-l ,MPI_C0MM_W0RLD ,&status) ; 

> 


> 

mpz_clear (N) ; 
mpz_clear (factor) ; 

f close (OUT) ; 

if (myrank==0)  //Compile  the  information  from  the  different  processors  into  a  single  file. 

{ 

lastOUT  =  f openC'MultipliersData. txt" , "w") ; 

for(i=0;  i<np;  i++) 

{ 

tens  =  i/10; 
f name [2]  =  ’O’+tens; 
f name  [3]  =  ’O’+i'/olO; 

IN  =  f open(fname , "r") ; 
single  =  fgetc(IN); 
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do 

{ 

fprintf  (lastOUT,  "#/0c"  ,  single)  ; 
single  =  fgetc(IN); 

}  while (single ! =E0F) ; 
fprintf (lastOUT, "\n") ; 
f close(IN) ; 
f flush (lastOUT) ; 

> 

f close (lastOUT) ; 

> 

MPI_Finalize() ; 
return  0; 


//Only  difference  between  this  and  the  standard  algorithm  (with  Fast  Return)  is  the  check  for  a  message 
//from  another  processor  indicating  completion, 
int  SQUFOF(mpz_t  N,  mpz_t  factor,  int  dDi  [3] ) 

{ 

const  int  STOP  [64]  =  {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31, 
32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63}; 
mpz_t  root ; 
struct  form  A,firstA; 
struct  form  saved [30] ; 
mpz_t  r ; 

int  found [NUMPROC] ; 

MPI_Request  waiting [NUMPROC] ; 
long  long  i ; 

int  pow=30,size=0,success=0,  stop=1000000000; 
int  np,myrank,f lag=0, j ,k; 

MPI_Status  status; 

MPI_Comm_rank(MPI_COMM_WORLD ,  &myrank) ; 

MPI_Comm_size (MPI_C0MM_W0RLD ,  &np) ; 

for(i  =  0;  i<np;  i++)  //Prepare  to  receive  messages  from  the  other  processes. 

{ 

if (i ! =myrank) 


MPI_Irecv (&f ound [i] , 1 , MPI.INT , i , STOP [i] , MPI_C0MM_W0RLD , &wait ing [i] ) ; 
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> 

//Initialize  memory  and  calculate  first  form. 
mpz_init (root) ; 
mpz_sqrt (root ,N) ; 

mpz_init_set_ui (A . QO , 1) ; 
mpz_init_set (A. P, root) ; 
mpz_init (A . Ql) ; 
mpz_init (A .b) ; 
mpz_init_set_ui (r , 1) ; 

reduce (&A, root ,N) ; 

mpz_init (f irstA . QO) ;  mpz_init (f irstA . P) ;  mpz_init (f irstA . Ql) ;  mpz_init (f irstA .b) ; 
ftof (&f irstA, A) ; 

cFracStep(&A,root) ; 
cFracStep(&A,root) ; 
i  =  0; 

//This  is  the  basic  step  in  the  algorithm. 

//Searches  the  given  segment  for  a  square.  Uses  Fast  Return  to  test  any  squares  found. 

//Checks  each  cycle  for  a  message  from  another  processor. 

while ( ( ( ! isSquare(A. Q0,r) ) I  I ( ! testF (A, saved, f irstA, i , r ,root ,N,&success ,f actor ,dDi) ) )&&(i<stop)&&( ! flag) ) 

{ 

if (i==pow) 

{ 

mpz_init (saved [size] . QO) ;  mpz_init (saved [size] . P) ; 
mpz_init (saved [size] .Ql) ;  mpz_init (saved [size] .b) ; 
ftof (&saved [size] ,A) ; 
size  =  size+1; 
pow  =  2*pow; 

> 

cFracStep(&A,root) ; 
cFracStep(&A,root) ; 

if ( ( !mpz_cmp_ui(A. QO, 1) ) I  I ( !mpz_cmp_ui (A . Ql , 1) ) ) 
i  =  stop; 
i  =  i+1 ; 

for(j  =  0;  ( j<np)&&( ! f lag) ;  j++) 
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{ 

if ( j ! =myrank) 

MPI_Test (&waiting[j] ,&f lag,&status) ; 
if (flag) 

{ 

if (dDi [1] ) 

printfC'Rank  0/0d  flag  by  70d\n"  ,myrank,  j)  ; 
for(k=0;  k<np;  k++) 

{ 

if ( (k ! = j ) && (k ! =myrank) ) 
MPI_Cancel(&waiting[k] ) ; 

> 

> 

> 

> 

//Cancel  unused  communications. 
for(j  =  0;  ( j<np)&&( ! f lag) ;  j++) 

{ 

if ( j ! =myrank) 

MPI_Test C&waiting [j] ,&f lag,&status) ; 
if (flag) 

{ 

for(k=0;  k<np;  k++) 

{ 

if ( (k ! = j ) && (k ! =myrank) ) 

MPI_Cancel(&waiting[k] ) ; 

> 

success  =  0; 

> 

> 

//Tell  the  other  processors  about  the  factor  found, 
if (success) 

{ 

for(k  =  0;  k<np;  k++) 

{ 

if (k ! =myrank) 

MPI_Cancel (&waiting  [k] ) ; 

> 
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if  CdDi  [1] ) 

printf  ("Found  by  rank  °/0d\n"  ,myrank)  ; 

> 

//Clear  memory  used, 
ford  =  0;  i<size;  i++) 

{ 

mpz_clear (saved [i] .Q0) ; 
mpz_clear (saved [i] .P) ; 
mpz_clear (saved  [i] .Ql) ; 
mpz_clear (saved [i] .b) ; 

> 

mpz_clear (A . Q0) ; 
mpz_clear(A.P) ; 
mpz_clear (A . Ql) ; 
mpz_clear (A.b) ; 

mpz_clear (f irstA . Q0) ; 
mpz_clear (f irstA.P) ; 
mpz_clear(f irstA. Ql) ; 
mpz_clear(f irstA. b) ; 

mpz_clear(r) ; 
return  success; 


> 


